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ABSTRACT 


Shape control of spaceborne antennas can provide the ability to correct for effects 
such as thermal distortion and manufacturing errors as well as control the shape of an 
antenna’s radiated beam. This thesis examines the performance of piezoceramic actuators 
in producing static deformation of a cantilever beam and analyzes the optimal placement 
of actuators to best approximate a desired deformation profile. Predictions of actuator 
effectiveness at producing curvature are developed using an Euler-Bernoulli model. An 
algorithm to determine the optimal locations and input voltages for a fixed set of actuators 
to achieve a desired deformation profile of a cantilever beam using embedded Nelder and 
Mead simplex search routines is presented and evaluated for two shape functions and 
various combinations of actuators. Experimental measurements show that the Euler- 
Bernoulli model provides a reasonable prediction of actuator performance at low input 
voltage but does not account for nonlinear behavior of the piezoceramic and the effects of 
hysteresis and transverse stresses. Further experiments demonstrate the ability of four 
piezoceramic actuators to produce an approximation of a parabolic deformation profile of 
a cantilever beam and illustrate the importance of considering these effects in determining 


the required actuator input voltages. 
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1. INTRODUCTION 


Á. BACKGROUND 


Control of the shape of long flexible structures is of interest for a variety of space- 
related applications. Many spacecraft designs include slender flexible appendages as 
components of on board antenna systems or sensors. Maintaining the desired shape of 
such an appendage can be a central requirement for obtaining optimal performance from 
its associated system. Effects such as thermal distortion, manufacturing errors, spacecraft 
acceleration and vibration can all cause deviation of the shape of a component from its 
desired profile. Active means to counter these and other effects are therefore of 
considerable interest. 

Piezoceramic actuators offer an attractive means to produce deformation of 
flexible structures. These actuators have no moving parts and require only a supplied 
electrical voltage to produce an electric field across the actuator material in order to 
function. Piezoceramics can be bonded to or, when properly insulated, embedded inside 
the structure they are intended to control. Applying an electric field to a piezoceramic 
actuator causes deformation of the piezoceramic material which, in turn, applies stress to 
the structure it 1s attached to. 

While substantial research effort has been devoted to the employment and 
placement of piezoceramic actuators for the active suppression of structural vibrations, 
considerably less attention has been focused on the use and optimal placement of 
piezoceramics for quasi-static shape control of flexible structures. Although 
piezoceramics are limited to the production of relatively small structural deformations, 
they may be more than adequate for certain applications, such as countering thermal 
effects or manufacturing errors or producing a small desired deformation profile in a 


flexible structure. 


B. SCOPE OF THESIS 


This thesis examines the performance of piezoceramic actuators in producing static 
deformation of a thin cantilever beam and analyzes the optimal placement of actuators to 
best approximate a desired deformation profile. Models of common configurations of 
piezoceramic actuators bonded to the surface of the beam are developed from the Euler- 
Bernoulli model of Crawley and Anderson [Ref. 1] to compare the effectiveness of the 
different configurations and examine the effects of actuator thickness, width and material 
properties on the ability of the actuator to produce curvature of the composite 
actuator/beam structure. Predictions of actuator performance from the Euler-Bernoulli 
model are compared to experimental measurements of the tip displacement and surface 
strain of a cantilever beam for several actuator configurations. 

An algorithm to determine the optimal locations and input voltages for a fixed set 
of actuators which minimize the error between the actual and desired shape of a cantilever 
beam under piezoceramic actuator induced deformation is presented. Optimization of the 
actuator locations and voltages is accomplished using the Nelder and Mead simplex search 
algorithm for the optimization of nonlinear functions. The MATLAB software package is 
used to implement and test the algorithm for two shape functions and various 
combinations of actuators. Capabilities and limitations of the algorithm are discussed. 
The predicted performance of four actuators placed in the optimum locations determined 
by the algorithm to approximate a parabolic beam deformation profile is compared to the 


results from a second set of experiments. 


II. PIEZOCERAMIC MATERIAL PROPERTIES 


Piezoceramics are part of a class of materials which exhibits piezoelectric behavior. 
Such materials have the property that material strain is produced in response to an electric 
field applied to the material in line with the its "poling" direction. Similarly, strain due to 
mechanical stresses on the material will produce an electric field in the material poling 
direction due to the accumulation of electrical charge on the material surfaces normal to 
the poling axis. This accumulated charge results in a voltage between the two surfaces. 
Figure 1 shows the standard coordinate system used to describe the orientation of a 
piezoelectric material. The poling axis is normally designated as the 3-axis, with the 1- 
and 2-axes indicating the principal material directions orthogonal to the poling axis. 

Piezoelectric materials are typically orthotropic, meaning that the material 
properties are uniform along each principal material axis, but the properties along one 
principal axis may not be the same as those along another principal axis. Piezoceramics, 
however, are generally transversely isotropic, meaning that the material properties are the 
same in all directions in the plane orthogonal to the poling axis, depicted by the 1-2 plane 


in Figure 1. 


Figure 1. Piezoceramic Material Coordinate System Conventions 


The strain behavior of piezoelectric materials in response to applied electric fields 
and the voltage produced by mechanical stresses makes piezoelectric materials suitable for 
use as both actuators and sensors in a variety of applications. The voltage between the 
two surfaces normal to the poling axis in both actuation and sensing applications requires 
that these two surfaces be electrically insulated from one another, hence piezoelectric 
materials are generally bonded to the surface of conductive metal structures with the 
material poling axis normal to the bond.  Piezoelectric materials can, however, be 


embedded in electrically insulating materials, including many composites. 


A. CONSTITUTIVE RELATION 


The axial strains in a piezoelectric material are functions of the normal stresses 
applied to the material and the electric field applied in the material poling direction. These 


strains are given by the material constitutive relation: 








ne = eee 
«E S sem 
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where the 3-direction is the poling direction of the piezoelectric material, the 1 and 2- 
directions are the principal material axes orthogonal to the poling direction, ¢ 1s the axial 
strain in the i-direction, Fi, Vj and ds; are the material elastic moduli, Poisson ratios and 
piezoelectric coefficients, respectively, o; is the axial stress applied in the i-direction and ¢3 
is the electric field applied in the material poling direction. The piezoelectric coefficients 
are generally dependent on the value of the applied electric field but are typically assigned 


a constant value which provides a linear approximation of the material's piezoelectric 


behavior over a specified operating range. The material elastic moduli are also field 
dependent and are typically specified for constant electric field (short circuit) and open 
circuit conditions along the poling axis. The moduli for constant electric field are used in 
Equation 2.1. Typical values of the material properties for various types of lead zirconate 


titanate (PZT), the most commonly used piezoceramic material, are listed in Table 1. 


PZT-4 PZT-5A PZT-8 PZT-5J PZT-5H 


E;, E; (x 10 N/m’) 
6.1 8.7 6.2 6.2 
6.9 9.9 7.1 7.1 
E3 (x 10'° N/m’) 
e Short Circuit 5.3 4.4 4.8 
e Open Circuit 12.6 10.6 n/a lil 


I 
8.2 
RD 
6.6 
ITAA 
verrons 


Table 1. Material Properties of Selected Piezoceramic Materials [After Ref. 2] 


















e Short Circuit 


• Open Circuit 











It should be noted that the piezoelectric constitutive relation is of the same form as 
the thermal-mechanical constitutive relation, with the piezoelectric coefficients 
corresponding to the material coefficients of thermal expansion and the electric field 
applied in the poling direction corresponding to the difference in material temperature 
from its reference temperature. This fact proves useful in modeling piezoelectric material 
behavior using finite element software which includes thermal-mechanical models but not 


piezoelectric models. 


B. ACTUATOR STRESS 


can be simplified to: 


If the piezoceramic material is under a state of plane stress in the plane orthogonal 
to the poling axis, i.e. оз = O, then the expression for the transverse strains in Equation 2.1 





(2.2) 


which can be solved for the transverse normal stresses: 
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where the reduced stiffnesses Qj are given by: 
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and the material Poisson ratios are related by: 


Р, 
Va si Vi 
1 


(2.5) 


Equations 2.1 through 2.5 apply to orthotropic piezoelectric materials and must be 
used for materials such as polyvinylidene flouride (PVDF) film which are not transversely 


isotropic in the plane orthogonal to the poling axis. Transversely isotropic materials such 


as PZT ceramics have the same properties in all directions in the 1-2 plane, i.e. 


Ë, = Е, 
Vi, = Va (2.6) 
1, = dz 


which allows Equation 2.3 to be further simplified to: 


Oif Qi Qo || £ da: 
Mi оГ wiy. 





(2.7) 
where: 
E 
(О = 2 
pp (2.8) 
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If the material is under a state of uniaxial stress along the 1-axis, i.e. o> = оз = О т 


Equation 2.1, then the expression for the axial strain along the 1-axis can be reduced to: 


О; 
сус О 
n 


(2.9) 


which can be solved for the normal stress along the 1-axis: 


санак) (2.10) 


The plot of Equation 2.10 in Figure 2 1s illustrative of the performance of a 
piezoelectric material as an actuator. An unconstrained piezoelectric material will have 
normal strain £; equal to the actuation strain d3;¢3 and will therefore produce no normal 
stress o;. In contrast, a piezoelectric material constrained to have no normal strain £; will 


produce the material's blocked normal stress: 


С ыо = Ез (2.11) 


In general, the normal stress in a piezoelectric actuator bonded to or embedded in a 
deformable structure will be somewhere between zero and the blocked normal stress, as 


determined by the material properties and geometry of the actuator and structure. 


O1 
-Е1дз1фз 


фз = constant 


d$; £i 


Figure 2. Axial Stress vs. Strain For Piezoelectric Material Under Uniaxial Stress 


HI EULER-BERNOULLI MODEL OF BEAM DEFORMATION DUE TO 
SURFACE-BONDED PIEZOCERAMIC ACTUATORS 


Piezoceramics can be used as actuators to produce curvature of a section of a 
beam. Voltage applied to a piezoceramic actuator bonded to the surface of a beam as 
shown in Figure 3 will induce normal strain of the same sign in both the 1- and 2- 
directions. The beam will resist expansion of the piezoceramic to its full actuation strain, 
producing normal stress in the piezoceramic which will cause curvature of the composite 
actuator/beam about the x- and z-axes. The curvature about the z-axis for a configuration 
such as that shown in Figure 3 will be significantly greater than the curvature about the x- 
axis due to the difference in the moments of inertia of the respective beam cross-sections. 
The resulting stress state in the actuator and the beam can be assumed to be plane stress 
provided that the material dimensions in the 3-direction are much smaller than the 
corresponding dimensions in the l- and 2-directions. The normal stresses in the 
piezoceramic for the plane stress condition are coupled in the 1- and 2-directions by the 
Poisson effect, as indicated by Equation 2.7. Solution of the coupled equations requires 
the use of a complex plate theory model, but a simple approximation of the curvature 


about the z-axis can be obtained using Euler-Bernoulli assumptions. 


3,y 


> 2 


d 


Figure 3. Piezoceramic Actuator Bonded to Beam 
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Euler-Bernoulli beam theory assumes that the transverse components of normal 
stress o» and o; are negligible compared to the axial component o; and hence uses 
Equation 2.10 to determine the normal stress in the piezoceramic along the beam axis. 
Comparison of Equations 2.7 and 2.10 reveals that the Equation 2.10 provides a lower 
estimate of the axial stress than Equation 2.7 for a given normal strain €; and applied 
electric field g3. The Euler-Bernoulli model therefore gives a conservative approximation 
of the bending moment produced by the actuator about the beam's z-axis for a given 
curvature of the beam about that axis. 

Crawley and Anderson [Ref. 1] compared two models of beam curvature due to 
induced strain actuation by piezoceramic actuators bonded to beam structures. The first 
model, initially presented by Crawley and de Luis [Ref. 2], assumed uniform axial strain in 
a pair of piezoceramic actuators bonded symmetrically to the outer surfaces of a beam 
undergoing actuator-induced bending. The second model assumed that the beam behaves 
as an Euler-Bernoulli beam, with a linear distribution of axial strain throughout the 
composite actuator-beam cross-section, as shown in Figure 4. The predicted beam 
curvature under both models was compared to that predicted by a detailed two- 
dimensional finite element model. The Euler-Bernoulli model was found to provide results 
that are within 0.1% of the finite element model over the full range of possible actuator- 
to-beam thickness ratios, while the uniform strain model results diverged significantly for 
actuator thicknesses greater than about 20% of that of the beam. The Euler-Bernoulli 
model was therefore judged to provide a much more accurate prediction of beam 
curvature produced by surface bonded piezoceramic actuators. 

Chaudry and Rogers [Ref. 4] performed a similar analysis of an actuator bonded to 
one side of a beam using a thermal-mechanical finite element model incorporating two- 
dimensional linear plane stress elements to model the effect of the actuator’s piezoelectric 
behavior. The Euler-Bernoulli and finite element models were again found to be in close 
agreement, differing by at most 3.9% for actuator-to-beam thickness ratios of 1 or less. 

Bronowicki and Betros [Ref. 5] refer to experimental work conducted by Griffin 


and Denoyer [Ref. 6] which demonstrated that the tip displacement of a beam with a pair 


10 


of piezoceramic actuators mounted symmetrically on opposite sides was 8% less than 
predicted by the Euler-Bernoulli model. The difference between the predicted and 
measured displacements was within the range of measurement and manufacturing 


tolerances of the experiment. 





Figure 4. Linear Strain Distribution for Euler-Bernoulli Beam 


All of the above references support the use of Euler-Bernoulli assumptions as a 
suitable model for analysis of the effectiveness of piezoceramic actuators in producing 
curvature of the composite actuator/beam. The development of this model by Crawley 
and Anderson [Ref. 1] was used as a basis for comparison of various actuator 


configurations in this thesis. 


A. GENERAL FORMULATION 


Crawley and Anderson [Ref. 1] expressed the total strain energy in an Euler- 


Bernoulli beam with attached or embedded piezoceramic actuators as: 


«Me es cen eb e е» 
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where / is the length, £ 15 the axial strain at y = 0, as shown in Figure 4, and x is the 


curvature of the actuator/beam section, which can be approximated by: 


к= (0 <<1) (372) 


The total area stiffness (EA)., the first moment of inertia (ES). and the second 


moment of inertia (E/), of the composite actuator/beam cross-section are respectively 


given by: 
(EA), = | E(y)w(y)dy (3.3) 
(ES), = | Ебђи(ујуду. (3.4) 
(EI), = ] Fowoy'& (3.5) 


where £(y) is the elastic modulus and w/(y) 1s the width of the composite actuator beam 
material as a function of the coordinate y. The first term inside the integral in Equation 
3.1, with the 1/2 multiplication factor, represents the strain energy per unit length for a 
composite beam with the specified cross-sectional properties composed of materials which 
exhibit linear stress-strain behavior in accordance with Hooke's Law. The second term 
and 1/2 multiplication factor represent a reduction in the strain energy per unit length due 
to the portion of the material strain which is induced by the piezoelectric effect rather than 
applied stress. 

The resultant force P4 and resultant moment M, due to the actuation strain A are 


respectively determined from: 


p? 


P, = |EWAWYOYHY (3.6) 


M, = | EQ)AQ)wO)yay (3.7) 


where the actuation strain A is determined by the actuator piezoelectric charge constant 
ds, the actuator thickness /, and the actuator input voltage V applied in the poling 


direction of the piezoceramic: 


imc 


3 


+ 
> (3.8) 
P 


If the actuator and beam properties are constant over the length of the 


actuator/beam section, Equation 3.1 can be integrated to give: 
1 
U = ;le (EA), +2e,x(ES), +? (EI), - 2P,£,-2M,x| (3.9) 


Equations for the axial strain £ at y=0 and the curvature < can then be obtained by 


differentiating Equation 3.9 with respect to & and « and using Castigliano's First Theorem 


E =Q, (i=1)2,...,n) (3.10) 


where q; is the ith displacement and Q; is the ith generalized external force. If no external 
forces or moments are applied to the beam, substitution of Equation 3.9 into Equation 


3.10 gives: 
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Ze "e. (E4), + 2x(ES), -2P,|=0 
(3.11) 


&, 
АЈ 1 
Es – 22, (Е5), +2«(EI),-2M,|=0 


which can be solved to obtain expressions for the axial strain at y = 0 and the curvature of 


the composite beam: 


Е P, (EI), P M, (ES), 
(БА) (КГ) - (ES); 

M, (ES) | (El), - M,(ES), — 
"7 (8), EAT ЈЕ e 





If the composite beam cross-section is symmetrical with respect to the z-axis, then 


the first moment of inertia of the cross-section is zero and Equations 3.12 simply to: 


E a 
” (EA) 
| G 19 
K EN, 
(ЕГ), 


The location of the neutral axis, which has zero axial strain as shown in Figure 4, 


can be determined by solving: 


ПЕ Е шк) О (3.14) 
for y, giving: 
EE (3.15) 
K 
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B. MODELS OF SELECTED ACTUATOR/BEAM CONFIGURATIONS 


Piezoceramic actuators are typically bonded to beams in one of four 
configurations: a single actuator, a pair of actuators bonded symmetrically to opposite 
sides of the beam, a stacked pair of actuators bonded to one side of the beam and stacked 
pairs of actuators bonded symmetrically to opposite sides of the beam. Prepackaged 
actuators, such as the PZT QuickPacks from Active Control eXperts, Inc., which have the 
piezoceramics and wiring encased in an insulating matrix material, can also be bonded to 
one or both sides of a beam. It is of interest to compare the effectiveness of each of these 
configurations in producing curvature of the composite beam. The general formulation of 
the Euler-Bernoulli model is applied to each of these configurations in the following 
subsections. | While the bond between an actuator and a beam is often negligibly thin 
when compared to the thicknesses of the beam and actuator, Baz and Poh [Ref. 7] noted 
that the thickness and material properties of the bond may contribute significantly to the 
composite beam’s mechanical behavior. The dimensions and modulus of the bond material 
are therefore included in the models for completeness. For simplicity, all piezoceramic 
actuators in the models with more than one actuator are assumed to have the same 


dimensions. 


1. Beam Deformation Due To Single Piezoceramic Actuator. 


Figure 5 shows a schematic representation of a piezoceramic actuator bonded to a 
section of a beam which extends in the +x-direction. The total area stiffness for the 
composite single actuator configuration is obtained by substitution of the elastic moduli E; 
and dimensions у; and /; into Equation 3.3, where no subscript, the subscript b and the 
subscript p denote the properties of the beam, bond and piezoceramic actuator, 


respectively: 
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(3.16) 





Side View Cross-Section 


Figure 5. Schematic of Single Piezoceramic Actuator Bonded to Beam 


The first and second moments of inertia of the composite beam cross-section are similarly 


obtained from Equations 3.4 and 3.5: 


| atts c 
(ES), = | Ewydy + | Ew, ydy + [Ew ydy 
- ; tb (3.17) 


pk і Ї 
= Ewa +4) + Ewé += 
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(EI), = вузу + [Ему з [Ewy dy 


= 
2 2 2 


TEN Deep 
- ESL) *Lh t (3.18) 





The resultant force P, and resultant moment M, due to an actuator input voltage V are 


obtained from Equations 3.6, 3.7 and 3.8: 


t 
gj btt, 
P= [E,Aw,dy = йүн (3.19) 
= 
7H, ; / 
М, == | E, Aw, ydy - Eds En y | (3.20) 


a, 
2 b 


Since the composite beam cross-section is not symmetrical with respect to the z- 
axis, an actuator input voltage will produce both non-zero axial strain £ at y = O and 
curvature x about the z-axis, as given by Equations 3.12, with the location of the neutral 


axis given by Equation 3.15. 


2. Beam Deformation Due To Symmetric Piezoceramic Actuator Pair 


Figure 6 shows a schematic representation of a symmetric piezoceramic actuator 


pair bonded to a section of a beam which extends in the -+x-direction. The total area 
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stiffness and first and second moments of inertia for a beam with a symmetric actuator pair 


can be computed as for the single actuator configuration using Equations 3.3 through 3.5: 


2E) E t zt e 
(ЕА) = [^ + fz, mdy + | Endy: [ndr Ts „Фу 
оя Jm е (3.21) 
2 2 2 2 2 


- Ewt * ZE,Ww,t, - 2E wt, 


Ф) а 
(ES), = [Ew ydy + А м, уду + jene 


[Iu ~ A 
2 2 E 
Lg, ott ote 
+ | E,w, ydy J Е w ydy (3.22) 
im 
2 2 
=0 


O 
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—|~+t, +t, 12. ) -— 
G Š 2 2 


(3.23) 
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Figure 6. Schematic of Symmetric Piezoceramic Actuator Pair Bonded to Beam 


The resultant force P4 and resultant moment M, as given by Equations 3.6, 3.7 and 3.8, 


depend on the input voltages V; to the actuator on the +y side and Vz to the actuator on 


the -y side of the beam: 


en is (2+) 
2 2 
Р, = [E,Aw,dy+ — [E,A,w,dy 
2+ [i 23! 
- E dw (V, *V,) 
LH +1, (5, | 
М, = [ЕА „уду + | E,A.w,ydy 
Dan _ я 
2 2 


Í і 
E dw, ha bs (V, V, ) 
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(3.24) 


(3.25) 


Since the composite actuator/beam cross-section is symmetric with respect to the 
z-axis, the axial strain at y = O and curvature about the z-axis are given by Equations 3.13. 


Substitution of Equations 3.24 and 3.25 into Equations 3.13 gives: 
ze E dw V, V, 
i (EA), 


t í 
Ар +1, + Hd -y,) 


(АІ) 


(3.26) 


к = 


c 


Thus, equal voltages applied to each actuator will produce axial strain with no curvature 
and input voltages of equal magnitude and opposite polarity will produce curvature with 
no axial strain at y — О. Input voltages of unequal magnitude will produce a combination 


of curvature and axial strain at y = 0. 


3. Beam Deformation Due To Stacked Pair of Piezoceramic Actuators 


Figure 7 shows a schematic representation of a stacked piezoceramic actuator pair 
bonded to a section of a beam which extends in the +x-direction. The total area stiffness 
and first and second moments of inertia for a beam with a stacked actuator pair are 
determined by substitution of the appropriate cross-sectional dimensions and properties 


into Equations 3.3 through 3.5: 


t 1 t t n 5 
2 z% PRA qu. =+21+21, 
(EA), = [ Ewdy+ [E,w,dy+ [E,w,dy+ [Ew,dy* — [E,w,dv 
А ш T t ie (3.27) 
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Side View Cross-Section 


Figure 7. Schematic of Stacked Piezoceramic Actuator Pair Bonded to Beam 
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The resultant force P, and resultant moment M, as given by Equations 3.6, 3.7 and 3.8, 


depend on the input voltages V; to the inner actuator and V; to the outer actuator: 


— ote 
E m dy + E Aw 
^ IP Wd „Је К „> (3.30) 
= E d w (V, +Ё,) 
Lagat, ть, 
= ЈЕ, Aw „ydy + ЈЕ, Aw ydy 
a ьн, (3.31) 


1 Д5 1 3t, 
= E d; W, 275 a Ko РА V, 


Since the composite actuator/beam cross-section is not symmetric with respect to 
the z-axis, the axial strain at y = O and curvature about the z-axis are given by Equations 
3.12. Actuator input voltages will, in general, produce both non-zero axial strain £ at y = 


O and curvature «x about the z-axis, with the location of the neutral axis given by Equation 


Э: 
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4. Beam Curvature Due To A Symmetric Pair of Stacked Piezoceramic 
Actuators 


Figure 8 shows a schematic representation of a symmetric pair of stacked 
piezoceramic actuators bonded to a section of a beam which extends in the +x-direction. 
The total area stiffness and first and second moments of inertia for a beam with a stacked 
actuator pair are determined by substitution of the appropriate cross-sectional dimensions 


and properties into Equations 3.3 through 3.5: 


t t t 
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Side View Cross-Section 


Figure 8. Schematic of Symmetric Pair of Stacked Piezoceramic Actuators Bonded to 


Beam 
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The resultant force P, and resultant moment M, as given by Equations 3.6, 3.7 and 3.8, 
depend on the input voltages V; and V2 to the inner actuators and V3; and V; to the outer 
actuators, where the odd numbered voltages are applied to the actuators on the +y side 


and the even numbered voltages are applied to the actuators on the -y side of the beam, 


respectively: 
[is +p) ey 
P,-  [E,Aw,dy«  [E,Aw,dy 
(Ls +2t, (ан, 
Stt + 2892, 
+ [EAw,dy* (Ели, (3.35) 
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(£52 uty] (E) 
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Since the composite actuator/beam cross-section is symmetric with respect to the 
z-axis, the axial strain at y = 0 and curvature about the z-axis are given by Equations 3.13. 


Substitution of Equations 3.35 and 3.36 into Equations 3.13 gives: 


— E da w (V, +V, +V, +V, ) 


° (FA) 


с 


Е d,,w Ë f, ele, у) о, Aye, An) y 
P P p 2 2 2 
(EI) 


к = 


c 


Hence, beam curvature with no axial strain at y = 0 can be obtained by setting V; = -V2 
and V;=-V,. Similarly, pure axial strain with no beam curvature is achieved when J; = |; 


and V3= V4. 


i Beam Deformation Due To Single ACX QuickPack Actuator 


Figure 9 shows a schematic representation of a single ACX QuickPack actuator 
bonded to a section of a beam which extends in the +x-direction. A QuickPack actuator 


consists of a stacked pair of piezoceramic wafers encased in a matrix material. 


26 










т T te 
SM z epa ииз uM dit 


Et 





Side View Cross-Section 


Figure 9. Schematic of Single ACX QuickPack Actuator Bonded to Beam 


The total area stiffness and first and second moments of inertia for a beam with a 
single ACX QuickPack actuator are determined by substitution of the appropriate cross- 


sectional dimensions and properties into Equations 3.3 through 3.5: 
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The resultant force P4 and resultant moment M, as given by Equations 3.6, 3.7 and 3.8, 


depend on the input voltages V; to the inner actuator and V; to the outer actuator: 
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Since the composite actuator/beam cross-section is not symmetric with respect to 
the z-axis, the axial strain at y = 0 and curvature about the z-axis are given by Equations 
3.12. Actuator input voltages will, in general, produce both non-zero axial strain £o at y = 


O and curvature «x about the z-axis, with the location of the neutral axis given by Equation 


ЗУБ. 


6. Beam Deformation Due To Symmetric Pair of ACX QuickPack 
Actuators 


Figure 9 shows a schematic representation of a symetric ACX QuickPack actuator 
pair bonded to a section of a beam which extends in the +x-direction. The total area 
stiffness and first and second moments of inertia for a beam with a symmetric ACX 
QuickPack actuator pair are determined by substitution of the appropriate cross-sectional 


dimensions and properties into Equations 3.3 through 3.5: 
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Figure 10. Schematic of Symmetric ACX QuickPack Actuator Pair Bonded to Beam 
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The resultant force P4 and resultant moment M, as given by Equations 3.6, 3.7 and 3.8, 
depend on the input voltages V; and V2 to the inner actuators and V; and V; to the outer 
actuators, where the odd numbered voltages are applied to the actuators on the +y side 


and the even numbered voltages are applied to the actuators on the -y side of the beam, 


respectively: 
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Since the composite actuator/beam cross-section is symmetric with respect to the 
Z-axis, the axial strain at y = 0 and curvature about the z-axis are given by Equations 3.13. 


Substitution of Equations 3.46 and 3.47 into Equations 3.13 gives: 
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Hence, beam curvature with no axial strain at y = 0 can be obtained by setting V, = -V; 
and V; — -V4. Similarly, pure axial strain with no beam curvature is achieved when V; = V; 


and V; = V4. 


c. ACTUATOR EFFECTIVENESS 


The results Euler-Bernoulli models of actuator produced curvature for the single, 
symmetric pair, single stacked and symmetric stacked pair configurations were compared 
for different thicknesses of an aluminum beam. All actuator configurations were assumed 
to cover the entire width of the beam and the bond between the piezoceramics and the 
beam was assumed to be negligibly thin. Material properties for Navy Type II PZT were 
used for the piezoceramic and a maximum applied electric field of 20 volts/mil, a typical 
operating limit, was applied to all actuators. Figure 11 shows the curvature of the 
composite actuator/beam as a function of the thickness of the individual piezoceramic 
elements for a 0.79 mm (0.031 in) thick beam. It can be clearly seen from Figure 11 that 
the symmetric actuator configurations are more effective at producing beam curvature 
than those with actuators on one side of the beam. This is due to the concentration of the 
piezoceramic at the outer edges of the composite beam for the symmetric configurations, 
where the actuator stress contributes most effectively to a bending moment about the 
central axis of the beam. An optimum piezoceramic thickness for each actuator type is 
also evident. Thicker than optimum actuators contribute more to the stiffness of the 
composite beam cross-section than to actuator bending moment and hence result in lower 
than optimum curvature. The optimum thicknesses for the two stacked configurations 
here are exactly half those of the corresponding single actuator configurations due to the 
assumption of negligible bond thickness. The stacked configurations therefore require half 
the input voltage of the single actuator configurations to achieve a given electric field 


across the piezoceramic to produce a desired level of actuation strain. 
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Figure 11. Beam Curvature vs. Piezoceramic Thickness and Actuator Type for 0.79 mm 
Aluminum Beam and 20 Volts/mil Electric Field applied to Actuators 


Figure 12 depicts the performance of the same actuators as in Figure 11 for an 
aluminum beam with twice the thickness. The optimum thickness of the actuators for each 
configuration 1s found to correspondingly double, however, the attainable curvature of the 
1.58 mm beam is only half that attainable for the 0.79 mm beam at the 20 volts/mil electric 
field limitation. This is due to the maximum actuation strain of the piezoceramic 
remaining constant. The limit on actuation strain for a given piezoceramic results in an 
inverse relationship between beam thickness and the curvature attainable using an actuator 
of optimum thickness. 

Figure 13 and Figure 14 show the same information as Figure 11 and Figure 12, 
respectively, but with the actuator thickness plotted as a percentage of the thickness of the 


aluminum beam. Figure 15 and Figure 16 show the effect of actuator width for a constant 
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actuator thickness of 0.26 mm and beam thicknesses of 0.79 mm and 1.58 mm 
respectively. Full width actuators are clearly the most effective in producing beam 
curvature for all actuator configurations. The reason for this can be seen by examining the 
curvature equation in Equations 3.26, or the equivalent equation for other actuator 
configurations. Actuator width appears as a factor of the entire numerator but of only one 
term in the denominator, given by Equation 3.23. Increasing actuator width therefore 


always results in an increase in beam curvature, provided all other terms remain constant. 
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Figure 12. Beam Curvature vs. Piezoceramic Thickness and Actuator Type for 1.58 mm 
Aluminum Beam and 20 Volts/mil Electric Field applied to Actuators 
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Examination of Equation 3.26 also reveals the material properties desirable in a 
piezoceramic actuator. High material modulus £,, piezoelectric coefficient d;, and 
operating voltage, or equivalently operating electric field, all result in increased actuator 


effectiveness. In general, the higher the product of these three properties, the more 


effective the actuator. 
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Figure 13. Beam Curvature vs. Piezoceramic Thickness as a Percentage of Beam 
Thickness and Actuator Type for 0.79 mm Aluminum Beam and 20 Volts/mil Electric 
Field applied to Actuators 
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Figure 14. Beam Curvature vs. Piezoceramic Thickness as a Percentage of Beam 
Thickness and Actuator Type for 1.58 mm Aluminum Beam and 20 Volts/mil Electric 
Field applied to Actuators 
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Beam Curvature at 20 Volts/mll 
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Figure 15. Beam Curvature vs. Piezoceramic Actuator Width and Type for 0.79 mm 
Aluminum Beam and 20 Volts/mil Electric Field applied to Actuators 
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Figure 16. Beam Curvature vs. Piezoceramic Actuator Width and Type for 1.58 mm 
Aluminum Beam and 20 Volts/mil Electric Field applied to Actuators 
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IV. ACTUATOR PLACEMENT 


A. DEFORMATION OF EULER-BERNOULLI BEAM UNDER ACTUATOR 
BENDING MOMENT 


| Transverse Deflection of a Beam Segment due to an Actuator Moment 


A section of laminated beam with uniform cross-sectional properties throughout its 
length will deform with uniform curvature over its length when subjected to a pure 
bending moment. For small displacements, defined as beam slope much less than 1, this 


behavior can be approximated by: 


dv, M, 
4^ (а), 





= ) (4.1) 


where x is the coordinate along the length of the beam, x; is the x coordinate of the start of 
beam segment k, /, is the x length of the beam segment, vo(x) is the transverse 
displacement of the beam's neutral axis in the y-direction, MM, is the bending moment 


applied to beam segment k about its neutral axis and 
(EN), = DEI, (4.2) 


is the flexural rigidity of the composite beam segment about its neutral axis. E, and /; are 
the Young's modulus and moment of inertia about the beam segment neutral axis of the ith 
lamina on segment А. 

The transverse displacement at any point x on beam segment А сап be obtained by 
integrating Equation 4.1 twice over the length of the segment. Integrating once gives the 


slope of the beam segment as a function of x: 
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0) = Sol we) ree ox) D < es x. +I, ) (4.3) 


Integrating a second time gives: 
M 
у, (х) = уо(х, ) + (6) 8) (x, SX A +.) (4.4) 
k 


Thus, for small displacements, the transverse displacement along the length of a beam 
segment is a function of the displacement and slope at the start of the beam segment, the 
curvature of the beam segment due to the actuator moment and the position x-x, on the 


beam segment. 


2: Transverse Deflection of Complete Beam due to Actuator Moments 
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Figure 17. Beam with Attached Piezoceramic Actuators 


Figure 17 shows a cantilever beam with n piezoceramic actuators bonded to the 


beam at locations x; (7 — 1,2,...,n). If the actuators are perfectly bonded to the beam over 
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their entire lengths and an actuator's bending moment is proportional to its input voltage, 


then each actuator will produce a beam curvature per unit voltage: 


M, к, 


Ke о 
m) V, 


(ЕМИ 07) (4.5) 


over its respective beam segment, where V, is the input voltage of actuator k and x; is the 
curvature of the corresponding composite actuator/beam section. A suitable value for К, 
for a given actuator/beam configuration can be obtained as discussed in Chapter II. 

If the deformation of the beam is static and no external loads are applied, then the 
deformation of the beam is strictly a function of the actuator locations, lengths, 
actuator/beam section curvatures per unit voltage and actuator input voltages. The static 
deformation of the beam on actuator k can be obtained by substituting Equation 4.5 into 


Equations 4.3 and 4.4, giving: 


VIE КОШЕ (у>, (x, SS +1,) | (4.6) 


v(x) -wn)en ex) (x (x, EX x, +1,) (АЛ) 


Similarly, since the sections of the beam between the actuators have no applied 
moments, and hence zero curvature, the static deformation of the beam between actuator 
k and К+ 1 can be obtained from Equations 4.3 and 4.4 with the moment terms set to zero. 


Thus, 
у, '(х) = vo (x, + 1,) (х, +1 <х< xn) (4.8) 
v (x) = v (x, + р) + Уа (20 + „х - (x, A) (x, nurse ху) (4.9) 
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Equations 4.6 through 4.9 include the terms vo(xj, vo'(xj), vo(xy* 1) and vo'(xy-- lj), 
which depend on the effects of actuators 1 through k-/ but not actuators k+/ through n. 
Thus, a complete expression for the transverse deflection of the beam as a function of x 
can be obtained using Equations 4.6 through 4.9 by starting at the root of the beam (x=0) 
and working outward to the tip (x=Z). 

Replacing the variable vo with y, the slope y' and displacement y of the beam are 


both zero prior to the first actuator: 


t dy 
c a) 
у(х) = 0 


The slope and displacement of the beam for x on the first actuator are functions of 


x and the properties of the first actuator only: 


y'(x) = yx) KV (x- x.) 
= 0-- KV (x — x,) 
= КУ(х-х) 


КУ, 
у(х) = y(x,)+ y'(x, (x -x,)4+ 2 ; (x- x) (x, < v x s) 








=0+0+ e (x-x,y 


AY 
i) 


The slope for x between the first and second actuators is the slope achieved at the 
end of the first actuator while the displacement for x in this range is the displacement at 
the end of the first actuator plus the slope at the end of the first actuator times the distance 


from the end of the first actuator: 
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y'(x) = ух, +h) 
= KV, 
у(х) x +1,)+y'(x, +/ 1-6 +1) 


saen (х, +2) 
| 
- KM x- x, > 


Continuing in this manner leads to the following general slope and transverse 


X, +Á <x<x) 





displacement relations for a beam with n attached piezoceramic actuators: 


a) For x prior to the first actuator 
'(x) = 0 
y e) (0€ x € x) (4.10) 
у(х) = 0 


b) For x on the first actuator 


y'(x) = KW, (x = x,) 








KV. осы dali 
у(х) = 1 (хх, ( 1 1 ) ( ) 
с) For x on actuator k 
k-1 
— ЕЕЕ ЕСЕ 
l (4.12) 
1 К = 2,3,...,п 
ж) 0) КИ s 3 


4) For x between actuator k and k+1 


(4.13) 


k 
у'(х) = A, E +1 <х< ны 


; 1 " И 
yx) B УКУ xx B^ A 
i=] 
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e) For x between actuator n and the free end of the beam (x = L) 


y'(x) = 23 КОР 
| 


z | (x, 1, < x< L) (4.14) 
у(х) = К | co i 
| 


B. OPTIMIZATION PROBLEM IN KUHN-TUCKER FORM 


It is desired to determine the optimum placement of piezoceramic actuators to 
approximate a desired beam shape function under actuator-induced deformation. 
Optimization requires a suitable cost function to serve as a measure of performance for a 
given actuator configuration. The optimum actuator configuration achieves the absolute 
minimum value of the cost function among all possible values of the input variables. For 
piezoceramic actuators, these input variables can be expressed as actuator position, length, 
input voltage, and beam curvature per unit actuator input voltage. 

An appropriate cost function for measuring actuator performance in static 


deflection of a cantilever beam is: 


Ј = fp)» (4.15) 


where x is the coordinate along the length of the beam and y(x) and y4(x) are the actual 
and desired transverse beam deflections at each point along the beam, respectively, as 
shown in Figure 17. Values of x range from 0 at the constrained end to L at the 
unconstrained end of the beam. Equation 4.15 is minimized by the combination of input 
variables that produces a y(x) which provides a continuous "least squares" fit to the 
desired shape function over the length of the beam. 


The actuator locations and lengths are constrained by: 


46 


x,20 

(От) 
ХЕ 

/>0 (ї=1,2,‚...‚п) 


(4.16) 


where 7 is the number of actuators and x; and /; are the locations and lengths of the 
actuators as shown in Figure 17. Equations 4.16 ensure that all actuators are located 
entirely on the beam with no overlap and that all have nonnegative length. Constraints on 
the maximum actuator lengths can also be included. Furthermore, the input voltages to 


the piezoceramic actuators must lie within an operating range specified by: 
ИИ Иоан) (4.17) 


where Vmn and Vina, are the actuator's minimum and maximum operating voltage and V; is 
the input voltage to the /th actuator. 
The problem of optimizing Equation 1 within the constraints of Equations 4.16 and 


4.17 can be expressed in Kuhn-Tucker (KT) form [Ref. 8] as: Minimize Equation 4.15 


subject to: 
в.(2) = –х, < 0 
g,(z)- x,, а a (i = 2,3,...,n) 
g_,(z)=x, +1, -Ls0 (4.18) 


ван (2) =-1, <0 (i =1,2,...,n) 
8 20+1+1 (2) = (У, = vA Xv, = Мы) < 0 (i = 2 A n) 


where the vector of input variables, z, is given by 
р И] (4.19) 
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A KT problem is typically solved by constructing a Lagrangian, given by: 
L- J(z) * и,8\(2) + и„ё.(г)+...+и„е„(т) (4.20) 


and solving the system of equations: 


a a), %0,,%0, ,, a.) 








P D а 9 

-0 ifg()«0 , o 
4, | ( = 12,...,m) 

>0 ifg,(z)=0 


where m is the number of constraint equations. Equations 4.21 with all 4 terms set to 0 
may have no solutions or one or more solutions for the values of the input variables, with 
each solution corresponding to a local minimum or maximum of Equation 4.15. A local 
minimum may or may not be the global minimum, if a global minimum exists. If Equation 
4.15 has either no finite global minimum or a global minimum outside the constrained 
range of input variables, the minimum value of Equation 1 over the allowable range of 
input variables may lie on a constraint boundary. The 4; terms in Equations 4.20 and 4.21 
allow for the possibility that the minimum value of the cost function in Equation 4.15 over 
the bounded range of possible input variables specified by Equation 4.18 lies on one or 
more constraint boundaries. The minimum value of the cost function on constraint 
boundary i can be found by including wu; as an additional variable to be solved for in the z 
vector and adding the corresponding constraint equation g;=0 to the system of equations 
to be solved in Equations 4.21. Differentiating Equation 4.20 with respect to the 


additional variable и; сап Бе seen to return the appropriate constraint equation. 
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C. SELECTED SHAPE FUNCTIONS FOR ANALYSIS 


An analytical expression for the cost function to be minimized can be obtained by 
substituting Equations 4.10 through 4.14 into Equation 4.15, provided that all terms 
involving the desired beam shape function y (Xx) are integrable. Two suitable shape 
functions were selected for evaluation of optimal actuator placement. The first, a 


parabolic shape function given by: 
y (v) = Cy (4.22) 
was selected as a due to its mathematical simplicity and its applicability to common 


structures such as antennas. This shape function has uniform curvature C over the length 


of the beam. A second shape function given by: 
y4G) = СП - соз(тх) | (4.23) 


where m = z/L, was chosen to analyze the effect of non-uniform desired curvature and an 


inflection point on the optimization process. This function has curvature: 


e" Сл? ЛХ 
ушык z s cod | (4.24) 


which is positive over the first half of the beam and negative over the second half. The 
inflection point for this function is at x = 0.5L, where the curvature is zero. 
Substituting Equation 4.22 and Equations 4.10 through 4.14 into Equation 4.15 


gives an expression for the cost function associated with the parabolic desired shape: 
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where the coefficients a; and b, are obtained by rewriting the beam displacements in 


Equations 4.10 through 4.14 as polynomials in x: 


For x prior to the first actuator: 


у(х) = 0 (0 ax =< x, } (4.26) 
For x on actuator X: 


x у tn 
(О a та ата, pa е ў (4.27) 
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For x between actuator k and K+ 1: 


Reet = x 
y(x) = bx + D E TA | ) (4.28) 
For x between actuator n and the free end of the beam (x = L): 
y(x) = b,x +,, (x, 17 = L) (4.29) 





giving: 
KY 
ay, = 2 | к= л) 
= (k — 1) 
" k-1 
TN - K,V,x, КИН (k = 2,3,...,n) 
1=1 
K V.x ° 
e i 
4 | (4.30) 


u K 
2 


k— 
> KV I, (x z (k = 205 
i=] 
k 
bye, (12 99) 
i=] 


k 
ВК (x x (k = 12; l) 


An equation for the cost function associated with the second shape function can 
similarly be obtained by substituting Equation 4.23 and 4.26 through 4.30 into Equation 
4.15, giving: 


5] 


2 


5 5 
(+) х, 
a 5 t 205,0, 


2 | MI 
J 2 C!| x, - Ssin(mx,) « 2 4 — sin(2mL 
E) 
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UNE 
+[2a,, (а, к тылы 
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+ 2а, (а, Ба oeta) z + (а, й б); I, 


^ а h I, )cos(m(x, +1, )- х5 cos(mx, )| 


I — I 2, sin(m(x, Lii ) B 2 u 5) E | 
А A t C [cos(m(x, +1, )) = cos(mx А ) 
+ "u^ [C A sin(m(x, d ) - x, sin(mx, | 


P Җы-С)с [sin(m(, +1,))—sin(me,) 
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b * 2b, (b, -C) 
+ (Бок е Ио — X, -1,) 
п-1 b C 
+ 2, i = |cos(mx, ..) = соз(т(х, = )| 
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(continued on next page) 
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4 TW P sin(mL) — (x, zi 1, )sin(m(x, +£, ) 
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”. ALGORITHM FOR SOLUTION OF OPTIMIZATION PROBLEM FOR 
ACTUATORS WITH FIXED DIMENSIONS AND PROPERTIES 


The MATLAB software package was used for development of an algorithm to 
optimize actuator placement and voltage for a given shape function and fixed actuator and 
beam dimensions and properties. The shape functions of Equations 4.22 and 4.23, with 
Equations 4.25 and 4.31, respectively as corresponding cost functions, were used for 
evaluation. Inspection of Equation 4.25, with appropriate substitutions from Equations 
4.30, reveals that its terms are up to fifth order polynomials in the actuator locations x; and 
lengths /; multiplied by second order polynomials in the actuator curvatures Kj; 
Equation 4.31 is similar, with sine and cosine terms replacing some of the polynomial 
terms. If the actuator dimensions and properties are fixed, then the lengths /; and 
curvatures per unit input voltage will be constant values. Either cost function is therefore 
minimized by optimizing the remaining variables x; and V;. Equivalently, if only the 
actuator lengths are specified, the cost function is minimized by optimizing the actuator 
locations x; and curvatures к; = К, The actuator configuration and cross-sectional 
properties can then be chosen to maximize the actuator curvature per unit voltage as 
discussed in Chapter III. 

The MATLAB Optimization Toolbox includes three functions intended for the 
optimization of multivariable functions. Two functions, fminu and fmins, perform 


unconstrained optimization of the input variables of a cost function, while one function, 
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constr, performs optimization of the input variables subject to constraints on their values. 
Fmins uses the simplex search algorithm developed by Nelder and Mead. Fminu uses the 
quasi-Newton method of Broyden, Fletcher, Goldfarb and Shanno (BFGS) with a mixed 
quadratic and cubic line search procedure to determine a search direction at each iteration. 
Constr also uses a BFGS quasi-Newton method to solve a quadratic subproblem at each 
iteration of a sequential quadratic programming routine. Grace [Ref. 9] provides a 
thorough description of the BFGS and sequential quadratic programming algorithms but 
omits detailed discussion of the Nelder and Mead technique. 

Initial attempts to optimize actuator placement and voltages simultaneously using 
constr and fminu proved unsuccessful. These algorithms returned unreliable results with 
minimal change in the initial actuator locations and voltages specified for each test, 
regardless of the settings of the tolerances on the input variables and cost function 
evaluations and other options associated with the functions. Fmins produced much better 
results in general but was found to return erroneous results for some sets of initial 
conditions for three or more actuators. The likelihood of erroneous results increased 
directly with the number of actuators evaluated. A more robust method of optimizing 
actuator locations was therefore judged to be needed. 

A revised approach to the optimization problem was developed, drawing on the 
work of Clark and Fuller [Ref. 10] and Wang, Burdisso and Fuller [Ref. 11] to optimize 
actuator locations for structural acoustic control. The latter group noted that the 
mathematical inconsistency in order of the voltages and actuator locations in their 
objective function formulation requires independent solution for the voltages and 
locations. Both groups therefore determined optimum actuator voltages for each set of 
actuator locations evaluated in their algorithms to determine optimum actuator placement. 
A similar approach can be applied to this problem by using an outer optimization 
algorithm to determine the actuator locations which minimize Equation 4.25 or 4.31 with 
an embedded algorithm to determine the optimum actuator input voltages for each set of 
actuator positions evaluated. The cost function for the outer algorithm thus becomes a 


fifth order polynomial in the actuator locations x; with sine and cosine terms if Equation 
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4.31 is used, while the cost function for the embedded voltage optimization algorithm 
becomes a second order polynomial in the actuator input voltages V,,. 

A two-stage solution algorithm to determine the optimum actuator locations and 
voltages for a given set of actuator and beam dimensions and properties was implemented 
as described above using MATLAB. The algorithm was tested using either constr, fminu 
or fmins to perform the required optimization. Both constr and fminu were again found to 
return unreliable results with almost no change in the specified initial conditions but the 
fmins simplex search routine was found to be highly reliable, converging to the same 
results for a given number of actuators and set of actuator and beam dimensions and 
properties, over a wide range of specified initial actuator positions. This is consistent with 
the observations of Parkinson and Hutchinson [Refs. 12, 13], who noted that the Nelder 
and Mead Simplex (NMS) algorithm has proven to be robust, although less efficient than 
some unconstrained optimization algorithms. A block diagram of the algorithm 
incorporating fmins is shown in Figure 18. The outer stage of the algorithm uses fmins to 
perform a simplex search to determine the actuator positions, starting from an initial 
estimate, which minimize the cost function to a set tolerance. The inner stage of the 
algorithm uses fmins to determine the optimum actuator input voltages for each set of 
actuator positions evaluated by the outer stage. The effect of actuator input voltage 15 
therefore completely transparent to the outer stage, while the actuator locations are fixed 
for each evaluation of the inner stage. 

The use of the unconstrained optimization algorithm fmins to solve the constrained 
actuator location and voltage optimization problems is acceptable provided that the 
optimum solution is at the global minimum of the cost function and does not violate any of 
the constraints. Since the integrand of Equation 4.15 1s always non-negative, the value of 
the cost function will also be non-negative and the cost function will have a global 
minimum. The algorithm presented here is therefore suitable for shape functions which do 
not result in actuator overlap and which require actuator input voltages within the 
allowable range. Complex desired shape functions with significant curvature variation 


over distances on the order of an actuator length or less might result in solutions with 
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actuator overlap, requiring that the NMS routine be modified to include constraints or 
replaced with another suitable constrained optimization algorithm for both voltage and 


position optimization. 


Select Initial Estimate of 
Optimum Actuator Positions 
Select Initial Estimate of 
Optimum Actuator Voltages 


Determine Optimum Actuator 
Locations using Nelder and 


Determine Optimum Voltages 
for Each Set of Positions 
Evaluated using Embedded 






Mead Simplex Algorithm Nelder and Mead Simplex 
Algorithm 
Output 
Results 


Figure 18. Actuator Placement Optimization Algorithm 


E. VALIDATION OF OPTIMIZATION ALGORITHM 


1. Results for Parabolic Shape Function 


The optimization algorithm was implemented using the MATLAB program 
actplace.m and its associated subroutines, which are presented in the appendix. The 
program was tasked with determining the optimum locations for between one and five 
actuators on a one meter long cantilever beam to best approximate a desired parabolic 
beam shape with a tip displacement equal to 0.1% of the beam length. The beam was 


specified to have the properties of 7075-T6 aluminum, a length of 1 meter, a width of 5.08 
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cm (2.00 in.) and a thickness of 1.575 mm (0.062 in.). The beam curvature per unit 
actuator input voltage was calculated for a symmetrically mounted pair of 0.26 mm (10 
mil) thick Navy Type II piezoceramic actuators covering 75% of the width of the beam. 
Each actuator length was fixed at 10% of the overall beam length. The algorithm was 
tested with at least four sets of arbitrarily selected initial actuator locations for each 
number of actuators to ensure convergence to the same optimum values. The initial 
configurations tested included all actuators clustered at the root of the beam, all actuators 
clustered at the center of the beam and actuators distributed over the length of the beam. 
An additional test was performed with the initial actuator locations within a few percent of 
beam length of the optimum locations found on previous trials for comparison of the 
number of floating point operations required to reach convergence. Initial actuator 
voltages were the optimum voltages for the initial actuator locations for each case, 
determined using fmins, except where noted. 

The results for each number of actuators are tabulated with respect to each set of 
initial conditions in Table 2 through Table 6. The optimization algorithm proved to have 
minimal sensitivity to the initial actuator locations specified for up to five actuators, with 
two exceptions denoted by asterisks in Table 5 and Table 6. Grouping the initial locations 
for four OS adjacent to one another at the center of the beam and using optimum 
voltages for this initial configuration as an initial estimate for all voltage evaluations 
prevented the inner voltage optimization algorithm from converging within the set 
maximum of 10,000 iterations at some point during the optimization process. Using the 
same actuator locations with initial voltage estimates of 0 produced erroneous results with 
some actuator overlap, as shown in the Run 2 column of Table 5. Similarly, grouping the 
initial locations for five actuators adjacent to one another at the root of the beam produced 
erroneous results with the first actuator overhanging the fixed end of the beam, as shown 
in the Run 1 column of Table 6. For both of these cases, slightly changing the initial 
conditions to allow 1 to 2% of beam length spacing between the initial actuator positions 
produced results consistent with those for the remaining sets of initial conditions tested, as 


shown in the Run 5 columns of the respective tables. The optimum actuator positions 
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obtained in the remaining cases for each number of actuators differed by at most 107, 
which was the set tolerance on the actuator locations for termination of the fmins function. 
The optimum deformed beam shape is plotted versus the desired parabolic shape for each 
number of actuators in Figure 19 through Figure 23. 

The results of the actuator voltage optimization algorithm were compared to the 
results of exhaustive searches of all possible combinations of discrete sets of actuator 
voltages, for the one, two and three actuator cases. The actuator voltages were varied 
discretely over specified ranges for actuator locations fixed at the optimum locations 
found by the placement optimization algorithm. Figure 24 through Figure 27 show plots 
of the variation of error with actuator input voltage. For all cases, the error due to the 
voltages obtained from the optimization was less than the smallest error obtained by the 
discrete search. Additionally, the voltages obtained from the optimization algorithm 
coincide with the locations of the minimum error in the corresponding plots. 

The results of the actuator placement optimization algorithms were similarly 
compared to the results of exhaustive searches of all possible combinations of discrete sets 
of actuator locations for the one, two and three actuator cases. The actuator locations 
were varied discretely over specified ranges and the optimum actuator input voltages were 
determined for each combination of actuator locations using the voltage optimization 
algorithm. Figure 28 through Figure 31 show plots of the variation of error with actuator 
location. As for the voltages, the error due to the actuator locations obtained from the 
optimization was less than the smallest error obtained by the discrete search and the 
locations obtained from the optimization algorithm coincide with the locations of the 
minimum error in the corresponding plots. Table 7 shows the results for three actuators 
which can be compared to the results of the optimal placement algorithm in Table 4. 

The comparison of the slopes of Figure 25 through Figure 27 and Figure 29 
through Figure 31 illustrates the relative importance of each of the actuators in 
determining the error between the actual and desired beam shapes. The location and input 
voltage, or equivalently curvature, of actuator 1, the actuator closest to the cantilevered 


end of the beam have the most significant effect. Similarly, the location and input voltage 
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of actuator 2 have a greater impact than those of actuator 3. This is consistent with the 
fact that each actuator affects the shape of the beam from its location to the beam's free 
end but has no effect on the shape of the beam between the cantilevered end and its 
location. Actuator 1 therefore affects the greatest portion of the beam's length, followed 


by actuator 2, then 3, etc. 


2. Results for Second Shape Function 


Following evaluation of the optimization algorithm's performance with the 
parabolic shape function, the algorithm was tasked with determining the optimum location 
for between one and four actuators to best approximate the beam shape given by Equation 
4.23 with a beam tip displacement again equal to 0.1% of the beam length. The beam and 
actuators were given the same dimensions and properties as in the previous tests. The 
algorithm was tested with at least four sets of initial actuator locations for each number of 
actuators to check convergence to the same optimum values. The initial configurations 
tested again included all actuators clustered at the root of the beam, all actuators clustered 
at the center of the beam, actuators distributed over the length of the beam and initial 
locations within a few percent of beam length of the optimum locations found on previous 
trials. Initial actuator voltages were the optimum voltages for the initial actuator locations 
for each case. 

The results of the tests for the shape function of Equation 4.23 are presented in 
Figure 32 through Figure 36 and Table 8 through Table 11. The optimization algorithm 
showed greater sensitivity to initial conditions for this shape function than for the 
parabola, particularly for the four actuator test, but again converged to essentially the 
same results for a variety of sets of initial conditions for each number of actuators. 
Several sets of initial conditions prevented the inner voltage optimization algorithm from 
converging within the set maximum of 10,000 iterations at some point during the 
optimization process. Grouping the initial actuator positions at the root of the beam 


resulted in the algorithm converging on local minima of the cost function for the two and 


39 


four actuator cases, as shown in the Run 1 column of the respective tables. While the 
result of the two actuator case for these initial conditions, with the first actuator 
overhanging the fixed end of the beam, is obviously erroneous, the result of the four 
actuator case is a valid local minimum. Comparison of the actuator locations which result 
in the global minimum of the cost function in Figure 35 and the actuator locations 
resulting in the local minimum in Figure 36 illustrates this. 

The curvature of the desired shape function, given by Equation 4.24, is 
decreasingly positive from the root to the midpoint of the beam and then increasingly 
negative from the midpoint to the free end of the beam. It can thus be qualitatively seen 
that actuators must be located in the regions of high desired curvature near the ends of the 
beam to most closely approximate the desired shape. Similarly, actuators in the low 
curvature region near the midpoint of the beam provide relatively little benefit in reducing 
the cost function. The combination of three actuators producing successively lower 
curvatures near the root and one actuator producing a high curvature in the opposite 
direction near the free end of the beam in Figure 36 provides a relatively good 
approximation of the shape function over its high-curvature portion near the root of the 
beam but a poorer approximation of the shape function over its high-curvature portion 
near the free end. In contrast, the combination of two actuators near the root and two 
actuators near the free end of the beam in Figure 35 provides a slightly poorer 
approximation of the shape function over its high-curvature portion near the root of the 
beam but a much better approximation of the shape function over its high-curvature 
portion near the free end, resulting in a lower overall value of the cost function than the 
arrangement in Figure 36. Grouping the initial actuator positions near the root of the 
beam prevented the optimization algorithm from moving more than one actuator past the 
inflection point in the desired shape function at x = 0.5L. However, the initial conditions 
for runs 4 through 6 in Table 11, which start two actuators on each side of the inflection 


point, allow the optimization algorithm to achieve the correct optimum positions. 
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3. General Observations 


The two shape functions considered here illustrate the importance of properly 
selecting initial conditions for the optimization algorithm. Allowing some space between 
initial actuator locations was proven to be beneficial, but not always required to obtain 
convergence. The presence of inflection points or significant variations in the curvature of 
the desired shape function can result in actuator locations producing a local minimum in 
the error cost function and must therefore be considered in selecting initial conditions. 
Knowledge of the approximate locations of the optimum actuator positions should be used 
if possible, and multiple sets of initial conditions should be evaluated to ensure the 
minimum obtained is the global minimum. Exhaustive searches of discrete sets of actuator 
locations can be helpful in selecting initial conditions, but such searches become 


increasingly computationally intensive as the number of actuators is increased. 
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Figure 19. Optimal Location for One 1096 Length Actuator for Parabolic Desired Shape 





Initial Actuator 0.2 0.25 
Location (m 













Voltage (V 
Location (m 
Voltage (V 





Table 2. Optimization Algorithm Results for One 1096 Length Actuator for Parabolic 
Desired Shape 
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Figure 20. Optimal Locations for Two 1096 Length Actuators for Parabolic Desired 
Shape 











Locations (m 0.1 05 0.65 0.5 

Voltages (V 61.55 -50.78 12.97 26.95 
Locations (m 0.5394 0.5394 0.5394 0.5394 
23.50 23.50 23.50 23.50 
Voltages (V) 24.85 24.85 24.85 24.85 


Table 3. Optimization Algorithm Results for Two 10% Length Actuators for Parabolic 
Desired Shape 
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Figure 21. Optimal Locations for Three 10% Length Actuators for Parabolic Desired 
Shape 
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(RunNumber [| 1 | 2 | 3 | 4 | 5 | 6 
Initial Actuator 0.0 0.35 0.20 0.1 0.05 0.075 
Locations (m) 0.1 0.45 0.45 0.4 0.35 0.370 
0.55 0.70 0.7 0.65 0.665 
Initial Actuator 31.45 19.40 14.69 16.89 
Voltages (V) 4.58 17:02 19.48 18.18 
19.22 17.94 18.51 























Optimum 
Actuator 
Locations (m 
Optimum 16.67 16.67 16.67 16.67 16.67 16.67 
Actuator 18.32 18.32 18.32 18.32 18.32 18.32 
Voltages (V 18.30 18.30 18.30 18.30 18.30 18.30 


10° 


Table 4. Optimization Algorithm Results for Three 1096 Length Actuators for Parabolic 
Desired Shape 








Initial Actuator 0.0 0.3 0.1 
Locations (m) 





0.05 
0.25 





0.28 
0.39 













Initial Actuator 
Voltages (V) 





Optimum Actuator 
Locations (m) 





Optimum Actuator 
Voltages (V) 







7.528 | 7.528 
Operations (10' 


Table 5. Optimization Algorithm Results for Four 10% Length Actuators for Parabolic 
Desired Shape 
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Figure 22. Optimal Locations for Four 1096 Length Actuators for Parabolic Desired 
Shape 
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Figure 23. Optimal Locations for Five 10% Length Actuators for Parabolic Desired 
Shape 
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Initial Actuator 0.0 

Locations (m) 0.1 
0.2 
0.3 
0.4 

Initial Actuator 

Voltages (V) 


Optimum Actuator 
Locations (m) 


Optimum Actuator 
Voltages (V) 


15.7 14.5 
Table 6. Optimization Algorithm Results for Five 1096 Length Actuators for Parabolic 
Desired Shape 
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Figure 24. Variation of Error with Input Voltage for One Actuator at Optimum Location 
for Parabolic Desired Shape 
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Error vs. Actuator Input Voltagos for 2 Actuators 


Run 1 


Actuator Locations: 0.122 m à 5394 m 





Voltage 1 {V} 


. Veltage 2 (V) 


Figure 25. Variation of Error with Input Voltages for Two Actuators at Optimum 
Locations for Parabolic Desired Shape 


7 Eror vs. Actuator 1 & 2 Voltage at Optimum Actuator 3 Voltage 


Run 1 


Actuator Locations: 0.0728 m 034678 m 0.8641 т 





30 25 
Voltage 1 (V) 


0 
Voltage 2 {V} 40 35 
Figure 26. Variation of Error with Actuator 1 and 2 Input Voltages for Three Actuators 


at Optimum Locations and Optimum Actuator 3 Input Voltage for Parabolic Desired 
Shape 
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x 10* Error vs. Actuator 2 & 3 Voltage at Optimum Actuator f Voltage 


Run 1 
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Figure 27. Variation of Error with Actuator 2 and 3 Input Voltages for Three Actuators 
at Optimum Locations and Optimum Actuator 1 Input Voltage for Parabolic Desired 
Shape 
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Figure 28. Variation of Error with Actuator Location for One Actuator for Parabolic 
Desired Shape 
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Error vs. Actuator Locations for 2 Actuators 
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Figure 29. Variation of Error with Actuator Location for Two Actuators for Parabolic 
Desired Shape 
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(RunNumer | 1 | 2 | 3 | 4 | 5 | 6 — 


Actuator 0.00- 
























Location 0.20 
Ranges (m) 0.25- 
0.45 
0.55- 















0. 15 


Spacing 0.05 0.025 0.0125 0.00625 0.005 
Between 

Discrete 

Actuator 

Locations (m 


























Discrete 0.07500 0.075 
Actuator 0.36875 0.370 
Locations for 0.66625 0.665 


Minimum Error 





Optimum 19.40 | 15.88 | 16.96 16.96 16.87 16.89 
Actuator 17.62 | 19.48 18.68 18.68 18.06 18.18 
17.94 | 18.50 | 18.13 18.13 18.28 18.25 


Table 7. Results of Discrete Actuator Placement Algorithm for Three Actuators for 
Parabolic Desired Shape 
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Error vs. Actuator 1 8 2 Location at Optimum Actuator 3 Location 
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Figure 30. Variation of Error with Actuator 1 and 2 Location for Three Actuators with 
Actuator 3 at Optimum Location for Parabolic Desired Shape 
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Error vs. Actuator 2 & 3 Location at Optimum Actuator 1 Location 


Run 2 


D Actuator Leeations: 004m 0.35 m 045m 
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Figure 31. Variation of Error with Actuator 2 and 3 Location for Three Actuators with 
Actuator 1 at Optimum Location for Parabolic Desired Shape 
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Figure 32. Optimal Location for One 1096 Length Actuator for Equation 4.23 Desired 
Shape 





Location (m 
Voltage (V 
Location (m 







Optimum Actuator 39.19 39.19 Eo S 39.19 
Voltage (V 


2.86 


Table 8. Optimization Algorithm Results for One 1096 Length Actuator for Equation 
4.23 Desired Shape 
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Figure 33. Optimal Locations for Two 10% Length Actuators for Equation 4.23 Desired 
Shape | 


Number 0] 1" | 2 | GNE NNI UNS 
0 
0 


1* 
Initial Actuator 0.0 4 0.25 0.1 0.1 
Locations (m 0.1 5 0.65 0.25 0.75 







Voltages (V 25.70 -190.42 -54.75 -12.17 -31.68 
Locations (m 0.0989 0.7441 0.7441 0.7441 0.7441 
Voltages (V 2121 -29.32 -29.31 -29.31 -29.30 


Table 9. Optimization Algorithm Results for Two 1096 Length Actuators for Equation 
4.23 Desired Shape 
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Figure 34. Optimal Locations for Three 10% Length Actuators for Equation 4.23 
Desired Shape 
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Initial Actuator 0.0 0.35 0.20 0.1 0.05 
Locations (m) 0.1 0.45 0.45 0.4 0.20 
0.2 0.55 0.70 0.7 0.70 





















Initial Actuator 
Voltages (V) 













Did Not 
Converge 


Optimum 


Actuator 
Locations (m 
Optimum 
Actuator 









Terror 10") —] 315 [3158 | | sis | 345 — 
[Operations (10) | 245 | 297 | | 151 | 084 - 


Table 10. Optimization Algorithm Results for Three 10% Length Actuators for Equation 
4.23 Desired Shape 


Run м Гб 


Initial Actuator 
Locations (m) 






















Initial Actuator 
Voltages (V) 





Optimum 
Actuator 
Locations (m) 





Optimum 
Actuator 
Voltages (V) 


(10° 


Table 11. Optimization Algorithm Results for Four 10% Length Actuators for Equation 
4.23 Desired Shape 
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10 Actual vs. Desired Beam Shapes 
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Figure 35. Optimal Locations for Four 10% Length Actuators for Equation 4.23 
Desired Shape 
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Figure 36. Actuator Locations Resulting in Local Minimum of Cost Function for 
Equation 4.23 Desired Shape 


F. EFFECT OF NUMBER OF ACTUATORS FOR A GIVEN TOTAL 
ACTUATOR LENGTH 


The optimization algorithm was run for a total length of actuators equal to 50% of 
the beam length divided into one, two, three, four or five actuators to analyze the effect of 
dividing a given total actuator length into increasing numbers of smaller actuators. The 
algorithm was tasked with determining the optimum location for the actuators on a 0.7875 
mm (0.031 in.) thick 7075-T6 aluminum cantilever beam to best approximate a desired 
parabolic beam shape with a tip displacement equal to 0.1% of the beam length. The 
beam curvature per unit actuator input voltage was calculated for a symmetrically 


mounted pair of 0.26 mm (10 mil) thick Navy Type II piezoceramic actuators covering 
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75% of the width of the beam. The results for each number of actuators are tabulated in 
Table 12 and plotted in Figure 37 through Figure 41. These results demonstrate that 
increasing the number of actuators making up the total actuator length contributes 
substantially to minimizing the error in approximating the desired shape. Doubling the 
number of actuators from one to two and from two to four each resulted in about an order 


of magnitude decrease in the error cost function. 


Actuators 
Length(s) (m 
Optimum 0.0852 0.0600 


Actuator 0.4684 
Location(s) (m) 


Optimum 
Actuator 
Voltage(s) (V) 





Table 12. Optimization Algorithm Results for 50% Total Actuator Length for Parabolic 
Desired Shape 
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Figure 37. Optimal Location for One 50% Length Actuator 
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Figure 39. Optimal Location for Three 16.67% Length Actuators 
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Figure 40. Optimal Location for Four 12.5% Length Actuators 
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Figure 41. Optimal Locations for Five 10% Length Actuators 
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V. EXPERIMENTAL ANALYSIS 


Experimental analysis was conducted in two parts. The first part consisted of 
comparison of the performance of several actuator configurations in producing beam 
curvature with the predictions from the corresponding Euler-Bernoulli models. The 
Euler-Bernoulli predictions neglected the effects of transverse stresses produced by the 
piezoceramic, assumed a linear strain distribution across the thickness of the composite 
actuator/beam and assumed a linear relationship between beam curvature and actuator 
input voltage. The second set of experiments examined the deformation of a beam due to 
four single actuators in the optimal locations to best approximate a parabolic displacement 


profile. 


A. COMPARISON OF ACTUATOR CONFIGURATIONS 


| Physical Configuration 


A schematic of the equipment setup for the experiment to compare the 
effectiveness of actuator configurations is shown in Figure 42.  Piezoceramic actuators 
were bonded to a 0.031 in. (0.79 mm) thick 7075-T6 aluminum beam which was 
cantilevered such that its length was horizontal and width was vertical to allow bending of 
the beam to take place in the horizontal plane. The single actuator, symmetric pair and 
stacked pair configurations shown in Figure 5, Figure 6 and Figure 7, respectively, were 
each tested on separate beams. Each beam had an overall length of 24 in., 2 in. of which 
was held in the clamp, and a width of 0.875 in. All actuators were 2.5 in. long, 0.75 in. 
wide and 0.26 mm thick Navy Type II piezoceramic. 

A XANALOG nonlinear dynamic systems modeling, simulation and real-time 
controls testing system was used to control actuator input and record experimental data. 
The XANALOG is a PC-based system incorporating an internal digital signal processor 


controlled by NL-SIM software running in a Windows environment. NL-SIM provides a 


E 


graphical interface which allows a control system model to be constructed using block 
diagrams and then compiled and operated in real time. System inputs and outputs can be 
displayed in both digital and graphical formats while the system is running and recorded to 


data files for post-experiment processing. 


Signal 
Conditioner 


Strain 
Gage Aluminum Beam 












Actuator 


Amplifier 
(15x) : 
| Laser 
Displacement 
Sensor 


Figure 42. Experimental Setup for Comparison of Actuator Configurations 


Beam tip displacement was measured at a distance of 20.125 in. from the center of 
each actuator using an NAiS ANL1651AC infrared laser analog displacement sensor with 
an output of 0.1 volts per millimeter of displacement. Beam surface strain was also 
measured directly opposite the center of the single and stacked pair actuator 
configurations in the beam length and width directions using Measurements Group Inc. 
CEA-13-125WT-350 bi-directional strain gages with a 2.16 gage factor. Each gage was 
connected to a Measurements Group Model 2120A Strain Gage Conditioner powered by 
a Model 2110A Power Supply. The Model 2120A output gains were set to 1732 to 


22. 


provide an output of 4 volts per 1000 microstrain measured by each strain gage. Analog 
to digital conversion of the laser and strain gage conditioner outputs was performed by an 
eight-channel CIO-DAS16/330 A/D board and digital to analog conversion of the actuator 
command signal from the XANALOG system was performed by a six-channel CIO- 
DDAO6 D/A board. Both boards were manufactured by Computer Boards, Inc. The +10 
volt analog output of the DDAOG was amplified by a 15:1 analog amplifier to provide a 
+150 volt range of input to each actuator. 

The command signal generator shown in Figure 43 was constructed using the NL- 
SIM software. The gain of 1.66667 was selected to allow integer command inputs to 
produce corresponding multiples of 25 volts input to the piezoceramic actuator when 
amplified by the analog amplifier. The second order filter shown was included to provide 
a gradual response to changes in the command input in order to prevent excitation of 
vibration in the beam. The filter's damping factor, С, was set to a value of 1 for critical 
damping and the filter’s natural frequency, @,, was set to 0.1047 rad/s, corresponding to a 
value of 1/60 Hz. The response of the filter allowed the XANALOG output to the analog 


amplifier to reach steady state to four significant figures approximately 1.5 minutes after a 
change in the commanded input was entered. 


Q ° Gain 
эшл 1.66667 
Set LEQ Se 


Figure 43. Block Diagram of XANALOG Actuator Command Signal Generator 
Implementation for Single Actuator Tests 
















Command 
Input 
(0, || 25 ,4,5,6) 


2. Experimental Procedure 


All of the experimental equipment was powered for at least 30 minutes prior to 
taking measurements to allow steady state conditions to be reached. The XANALOG 


system was configured to sample commanded actuator input and laser displacement sensor 
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and strain gage output at one second intervals to allow recorded data to be averaged in 
order to minimize the effects of noise in the sensor outputs. The laser and strain gage 
conditioner outputs were adjusted to indicate an average zero value as measured by the 
XANALOG system prior to the start of each measurement run. Each measurement run 
consisted of recording sensor outputs at actuator inputs of 0, 25, 50, 75, 100, 125, 150 
and again 0 volts. Each actuator input was allowed to remain steady for 60 seconds to 
provide 60 sensor output data points for each actuator input to be averaged. The next 
voltage level was commanded immediately after each 60 second measurement and allowed 
two minutes to reach steady state before the next 60 second measurement was initiated. 
At least three measurement runs were conducted for each experiment to verify 
repeatability of the results. Data for each measurement run was saved to a MATLAB 
.MAT data file for processing and plotting. Measurements were taken for the single 
actuator acting in both compression and tension to compare the effectiveness of 
unsymmetric actuators for these two cases. Measurements for the stacked actuator pair 
were taken for the actuator acting in compression, while measurements for the symmetric 


actuator pair were taken with one actuator acting in tension and the other in compression. 


3 Experimental Results 


Measurement data for the single actuator beam is presented in Table 13 and 
plotted in Figure 44 for the actuator acting in compression and Table 14 and Figure 45 for 
the actuator acting in tension. Data for the symmetric actuator and stacked actuator pair 
experiments are presented in Table 15 and Figure 46 and Table 16 and Figure 47, 
respectively. Positive strain indicates tensile strain for the single actuator acting in tension 
and compressive strain for all other experiments. All runs indicated in the tables and 
figures were conducted consecutively with less than five minutes between the completion 
of one run and the start of the next. The first run for all experiments except the single 
actuator acting in compression was the first use of that actuator on the day of the 


experiment. The first measurement run for the single actuator acting in compression was 
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conducted immediately after preliminary testing of the experimental equipment which used 
the same actuator. 

Significant hysteresis, or residual strain after the removal of applied voltage from 
the actuator, was apparent on the first run for all but the single actuator compression 
experiment. Udd [Ref. 14, p. 503] notes that hysteresis occurs in piezoceramic materials 
as the result of "energy dissipation due to internal sliding events in the polycrystaline 
piezoelectric body” and typically ranges from 0.1 to 10% of the strain due to an applied 
electric field prior to the removal of that field. Residual beam tip displacements of 31%, 
23% and 36% of the displacement at maximum input voltage were noted here for the 
single, symmetric and stacked actuator beams respectively. The hysteresis effect is much 
less apparent in the following measurement runs for each actuator since the zero- 
displacement reference was adjusted at the start of each run. Measurements for the single 
actuator acting in compression and runs 2 through 4 of the other actuator experiments are 
therefore displacements and strains above any hysteresis level in the material at the start of 
the run. The displacement measurements for these runs demonstrated excellent 
repeatability for all of the actuators. The repeatability of strain measurements for all but 
the single actuator beam was also excellent. Strain measurements for the single actuator 
showed somewhat more variation from run to run but still reflected the same general trend 
as the other measurements. 

The dashed lines in Figure 44 through Figure 47 depict the linear variation of beam 
tip displacement and beam surface strain predicted by the Euler-Bernoulli model, given by 


Equation 4.14 for one actuator: 


l 
Убе) КУ - X = (x, +Í 901, (5.1) 


l 
where Ë X, ==) is the distance between the center of the actuator and the 


measurement point at the beam tip. It can be seen from the figures that the Euler- 
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Bernoulli model provides a good approximation of the actual behavior of each 
actuator/beam combination at low actuator input voltage but differs increasingly from the 
actual behavior as voltage is increased. The actual relationship between beam tip 
displacement and actuator input voltage can be seen to be nonlinear. This is consistent 
with the fact that the published piezoelectric coefficient, d;;, of the piezoceramic 
represents a linear approximation of the material behavior over a specified operating 
range. The dotted lines in Figure 44 through Figure 47 represent second-order least 


squares curve fits of the form: 
(И) = сИ? + сИ + с, (5.2) 


to the observed displacements as a function of voltage. Coefficients determined for these 
curve fits and predicted values of c, from Equation 5.1 are given in Table 17. It should be 
noted that the values obtained for c? are negligibly small in magnitude. An expression for 
a second-order curve fit of actuator curvature as a function of actuator input voltage can 


be obtained by dividing Equation 5.2 by the actuator length /; and the measurement 


l 
distance Е cu in and omitting the negligible constant term: 


с +c V 


5-5 
| X — X, 73 


Values obtained for the ko and k; coefficients and the predicted curvature per unit voltage 


K(V) = -ky) «kV (5.3) 


are given in Table 18. 

The percentage of actual longitudinal strain in relation to the predicted strain was 
roughly 80% of the percentage of actual beam tip displacement in relation to the predicted 
displacement for the unsymmetric actuator configurations and roughly 140% for the 


symmetric actuator configuration. These percentages are too great to be explained by 
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measurement inaccuracies and variation of the beam and actuator material dimensions, 
which were determined to be within 2% of nominal values. Since strain was measured at 
the surface of the aluminum beam for the unsymmetric configurations and at the surface of 
the piezoceramic for the symmetric configuration, these results suggest that the strain 
distribution was not truly linear through the composite beam cross-section as shown in 
Figure 4. Rather, the strain at the aluminum beam surface was less than predicted by the 
assumed linear strain distribution but greater than predicted at the piezoceramic surface. 
The beam tip displacement measurements nonetheless provide a measure of the overall 
effect of the strain distribution over the entire actuator section of the beam at producing 
curvature of the beam’s neutral axis. 

Measured transverse strain tended to exceed the measured longitudinal strain by 
between 10% and 20%. The longitudinal strain had been expected to be significantly 
greater than the transverse strain due to the significantly larger moment of inertia of the 
beam cross-section about its longitudinal axis than about its transverse axis. The larger 
longitudinal moment of inertia should result in significantly less longitudinal curvature, and 
hence transverse strain, produced by the actuator. The strain gage conditioners for the 
longitudinal and transverse strain gages were swapped to eliminate the possibility of 
measurement equipment contributing to this result but found not to be a causing factor. 
The magnitude of the transverse strain measurements indicate that significant curvature 
about the beam’s longitudinal axis does in fact occur at least locally near the center of the 


actuator. 
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25 V 0.57(112) | 0.57(112) | 0.57 (113) 
50 V 1.20(118) | 1.18(116) | 1.19 (116) 
75 N 1.92 (125) | 189 (124) | 1.90 (124) 
100 V 2.65(130) | 2.65(130) | 2.65 (130) 
125 V 3.47(134) | 3.47(134) | 3.46 (134) 
150 V 4.38 (143) | 437(143) | 4.32 (141) 
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Table 13. Beam Tip Displacement and Strain Measurements for Single Actuator Acting 
in Compression 
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Figure 44. Beam Tip Displacement and Longitudinal Strain Measurements Single 
Actuator Acting in Compression 
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Table 14. Beam Tip Displacement and Strain Measurements for Single Actuator Acting 
in Tension 
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Figure 45. Beam Tip Displacement and Longitudinal Strain Measurements for Single 
Actuator Acting in Tension 
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Table 15. Beam Tip Displacement and Strain Measurements for Symmetric Actuator Pair 
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Figure 46. Beam Tip Displacement and Longitudinal Strain Measurements for 
Symmetric Actuator Pair 
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Table 16. Beam Tip Displacement and Strain Measurements for Stacked Actuator Pair 
Ácting in Compression 
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Figure 47. Beam Tip Displacement and Longitudinal Strain Measurements for Stacked 
Actuator Pair Acting in Compression 
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Table 17. Coefficients of Curve Fits for Observed Beam Tip Displacement as a Function 


of Actuator Input Voltage 
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Table 18. Coefficients of Curve Fits for Observed Actuator Curvature as a Function of 
Actuator Input Voltage 


B. BEAM SHAPE CONTROL WITH FOUR ACTUATORS 


1. Physical Configuration 


A schematic of the equipment setup for the four-actuator shape control experiment 
to compare the effectiveness of actuator configurations is shown in Figure 48. Eight 
piezoceramic actuator patches were bonded in four groups of two to one side of a 0.031 
in. (0.79 mm) thick 7075-T6 aluminum beam which was cantilevered such that its length 


was horizontal and width was vertical to allow bending of the beam to take place in the 
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horizontal plane. The beam was also supported at approximately two-thirds of its length 
by an air pad riding on a granite table. The beam had an overall length of 45 in., 2 in. of 
which was held in the clamp, and a width of 1.625 in. All actuator patches were 2.5 in. 
long, 1.50 in wide and 0.26 mm thick Navy Type II piezoceramic, hence each of the four 
actuator groups had a total length of 5 in. The actuators were placed at locations on the 
beam shown in Table 19 which were determined by the optimization algorithm to best 
approximate a parabolic deflected beam shape with four 5 in. long actuators on a 43 in. 
long beam. The XANALOG system was again used to control actuator input and record 
experimental data and the ANL1651AC laser displacement sensor was used to measure 
displacement at selected points along the length of the beam. Beam surface strain was 
also measured directly opposite the center of the actuator patch nearer to the cantilevered 
end of the beam in each of the four actuator groups using a Measurements Group Inc. 
CEA-13-125WT-350 bi-directional strain gage opposite the first actuator and a CEA-13- 
125UW-350 unidirectional strain gage opposite the other three actuators. Strain gage 
output was obtained with a Gould Model 56-1301-00 DC/Bridge/Transducer Signal 
Conditioner for each gage. DAS16/330 and DDAO6 boards were used for A/D and D/A 


conversion, respectively. Each +10 volt analog output channel of the DDAO6 was 


actuator. 










Distance from Percentage of Actuator 3 








Cantilevered End Curvature for Parabolic 


(in/cm) 





amplified by a 15:1 analog amplifier to provide a +150 volt range of input to each 
Shape Approximation 


11.52/29.26 100.0 








21.43/54.43 100.0 


Jable 19. Actuator Locations and Curvature Ratios for Parabolic Shape Approximation 
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Figure 48. Setup of Four Actuator Shape Control Experiment 


23 Experimental Procedure 


A series of measurements was initially taken to characterize the relationship 
between actuator curvature and input voltage for the single actuator configuration on the 
beam used for this set of experiments. Beam tip displacement measurements were taken 
with voltage applied to the first actuator only and then with voltage applied to the second 
actuator only using the procedures described for the previous set of experiments. The 
beam tip displacement measurement point was 38 in. from the center of the first actuator 
and 28.13 in. from the center of the second actuator. 

Two shape control experiments were then performed. Each experiment consisted 
of applying voltages predicted to produce the optimum approximation of a parabolic 


deformation profile along the beam’s length from its initial undeformed shape. Voltages 
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for each experiment were calculated to produce curvature of the four actuators in the 
proportions listed in Table 19. Identical voltages were applied to actuators 2 through 4 
since there was negligible difference between the curvatures required from these actuators. 
Corresponding input voltages for actuator 1 were determined for the first experiment by 
assuming a linear relationship between curvature and voltage as developed in Chapter II. 
Actuator 1 voltages for the second experiment were determined by solving Equation 5.3 
for the curvature produced by each actuator 2-4 input voltage and then solving Equation 
5.3 for the voltage to produce actuator 1 curvature equal to 91.4% of the actuator 2-4 
curvature. Coefficients kp and k; for Equation 5.3 of 5.19x10*? rad/m/V? and 6.2610" 
rad/m/V, respectively, were obtained from a curve fit of beam tip displacement to actuator 
l input voltage as described in the following Experimental Results section. 

Beam displacement was measured using the laser displacement sensor at locations 
of 18 in., 30 in. and 42.22 in. from the cantilevered end of the beam for each of the two 
experiments. It should be noted that the laser measurement axis was oriented 
perpendicular to the longitudinal axis of the undeformed beam. Increasing beam 
displacement therefore resulted in the displacement measurement points slightly further 
down the length of the beam than the nominal measurement locations. However, the small 
magnitude of the measured displacements relative to the length of the beam for these 
experiments made this effect insignificant. 

Three sets of measurements were taken for each of the two experiments at each 
measurement point. Each measurement run consisted of recording sensor outputs at 
actuator inputs given for the input stages in Table 20 and Table 21 for the first and second 
experiment, respectively. Each actuator input was allowed to remain steady for 60 
seconds to provide 60 sensor output data points for each actuator input to be averaged. 
The next input voltage stage was commanded immediately after each 60 second 
measurement and allowed two minutes to reach steady state before the next 60 second 
measurement was initiated. Data for each measurement run was saved to a MATLAB 
MAT data file for processing and plotting. All measurements were taken with the 


actuators acting in compression. The actuators were cycled through the input voltages 


109 


corresponding to one measurement run prior to the recording of data to minimize the 


effects of hysteresis on the repeatability of the measurements. 


Actuator 2,3,4 


Input Voltage (V) 


Actuator 2,3,4 


Input Voltage (V) 





Table 21. Nonlinear Model Actuator Input Voltages for Shape Control Experiment 2 


3. Experimental Results 


Displacement and stram measurement data for beam deformation due to the first 


actuator is presented in Table 22 and plotted in Figure 50. Data for beam deformation due 
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to the second actuator is presented Table 23 and Figure 52. Positive strain values indicate 
compressive strain. All runs indicated in the tables and figures were conducted 
consecutively with less than five minutes between the completion of one run and the start 
ofthe next. The first run for each actuator was the first use of that actuator on the day of 
the experiment. Some hysteresis was observed on all measurement runs but most apparent 
on the first run for each actuator, with residual beam tip displacements of 16% and 24% of 
the displacement at maximum input voltage occurring on the first run for the first and 
second actuator, respectively. Repeatability of the displacement measurements for both 
actuators and strain measurements for the first actuator was good for runs 2 through 4, 
but the strain measurements for the second actuator were not as reliable. Strain 
measurements for runs 1] and 3 of the second actuator were substantially lower than 
corresponding measurements for other runs of the first and second actuator. The final 
zero input strain measurements for both of these runs indicated residual tensile rather than 
compressive strain, as denoted by the asterisks in Table 23. The combination of low 
measurement values and final tensile strain measurements suggest that either a shift in the 
strain gage conditioner zero strain output voltage occurred shortly after the start of each 
of these runs or the bond between the actuator 2 strain gage and the beam may have been 
poor, leading to inconsistent results. 

The dashed lines in Figure 50 and Figure 52 represent the relationship between 
actuator input voltage and beam tip displacement predicted by the linear Euler-Bernoulli 
model, while the dotted lines represent second order least squares curve fits of the data 
from runs 2 through 4 for each actuator. Coefficients of Equation 5.2 for each of the 
curves are listed in Table 24 and the corresponding Equation 5.3 coefficients obtained are 
listed in Table 25. The ko and k; coefficients obtained for the two actuators agreed to 
within 6% and 2%, respectively. 

It can be seen from the figures that the Euler-Bernoulli model again provides a 
reasonable approximation of the behavior of the beam as a function of input voltage, as 
was observed for the previous beams tested. However, unlike the previous experiments, 


the Euler-Bernoulli model overpredicted the tip displacement of the 1.625 in. wide beam 
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for actuator input voltages below about 75 volts and underpredicted the tip displacements 
at higher voltages. Comparison of the percentages of predicted displacement obtained for 
the 1.625 in. wide beam in Table 22 and Table 23 with the corresponding percentages of 
predicted displacement obtained for the 0.875 in. wide beam in Table 13 and Table 14 
reveals that the percentage of predicted displacement obtained was always less for the 
wider beam. The reason for this is that the actuator produces curvature of the beam about 
both its transverse and longitudinal axes, a fact supported by the significant transverse 
strain measurements obtained. Curvature about the longitudinal axis increases the 
moment of inertia of the beam's cross-section, as shown in Figure 51, resulting in a 
greater resistance of the beam to bending about its transverse axis. The increase in 
moment of inertia is more significant for the wider beam than for the narrower beam. 

The results of shape control experiment 1 are tabulated in Table 26 and plotted in 
Figure 52 and Figure 53. The “linear model predicted” column of Table 26 and the 
dashed lines on the figures represent the desired displacement of the beam along its length 
to approximate a parabolic shape predicted by the Euler-Bernoulli model of actuator 
curvature and Equations 4.10 through 4.14. The “nonlinear model predicted” column in 
the table and dotted lines on the figures represent the displacement of the beam predicted 
by Equations 4.10 through 4.14 with the curvature terms K;V; replaced by Equation 5.3 
with the empirical coefficients kọ and k; equal to those obtained from the experiment to 
characterize the first actuator given in Table 25. The measured displacements displayed 
excellent repeatability and can be seen to agree much more closely with the beam shape 
profile predicted using the nonlinear empirical model of actuator curvature than for the 
profile predicted by the linear model. As for the single actuator tests, the linear model 
overestimates the actuator performance for the low input voltages and underestimates 
actuator performance at higher voltages. 

The results of shape control experiment 2 are presented in Table 27, Figure 54 and 
Figure 55. The "nonlinear model predicted" column in the table and dotted lines on the 
figures represent the desired displacement of the beam to approximate a parabolic shape 


determined using Equations 4.10 through 4.14 with the curvature terms K;V; replaced by 
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Equation 5.3, with the empirical coefficients kọ and k; equal to those for the first actuator 
in Table 25. The "linear model predicted" column in the table and dashed lines on the 
figures represent the displacement of the beam predicted by the linear Euler-Bernoulli 
model of actuator curvature and Equations 4.10 through 4.14 for the specified input 
voltages. The measured displacements at the 18 in. and 30 in. measurement points 
displayed excellent repeatability, but the displacements at the 42.22 in. measurement point 
showed somewhat more variation. The latter measurements were the first series taken for 
this experiment and therefore the most affected by hysteresis. Residual displacements of 
about 6% of the displacement at maximum input voltage were observed for runs 1 and 2 
at the 42.22 in. measurement point, but less than 2% residual displacement was observed 
on run 3 and for the measurements at the other measurement points. With the exception 
of runs 1 and 2 at the 42.22 in. measurement point, the measured displacements displayed 
much better agreement with the desired beam shape profile predicted using the nonlinear 
empirical model of actuator curvature than with the profile predicted by the linear model. 
The results of the shape control experiments demonstrate the importance of 
considering the nonlinear behavior of the piezoceramic in determining the appropriate 
actuator inputs for precision shape control applications. The significant effects of 
hysteresis observed indicate that this effect must also be considered. Actual shape control 
applications should incorporate some form of feedback of the beam’s shape to ensure that 


the desired deformation profile is attained. 
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Table 22. Beam Tip Displacement and Strain Measurements for First Actuator Acting in 
Compression 
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Figure 49. Beam Tip Displacement and Longitudinal Strain Measurements for First 
Actuator Acting in Compression 
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Table 23. Beam Tip Displacement and Strain Measurements for Second Actuator Acting 
in Compression 
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Figure 50. Beam Tip Displacement and Longitudinal Strain Measurements for Second 


Actuator Acting in Compression 
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Table 24. Coefficients of Curve Fits for Observed Beam Tip Displacement as a Function 
of Actuator Input Voltage 
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Table 25. Coefficients of Curve Fits for Observed Actuator Curvature as a Function of 
Actuator Input Voltage 



























Figure 51. Graphic Depiction of Deformation of Beam Cross-Section due to 
Piezoceramic Actuator Acting in Compression 
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Table 26. Beam Displacement Measurements for Shape Control Experiment I (Linear 
Model Voltages) 
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Figure 52. Beam Displacement Measurements for Shape Control Experiment 1 (Linear 
Model Voltages) (1,2 of 3) 
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Figure 53. Beam Displacement Measurements for Shape Control Experiment 1 (Linear 


Model Voltages) (3 of 3) 
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Beam Displacement at Linear Nonlmear 
Measurement Locations - mm Model Model 
Predicted | Predicted 

46.5/50/50/50 Volts Input 

х = 18.0 in 

x — 30.0 in 

x = 42.2 in 
93 .6/100/100/100 Volts Input 

x = 18.0 in 

x = 30.0 in 

x = 42.2 in 
141.1/150/150/150 


Volts Input 
х = 18.0 m 
x = 30.0 in 
x = 42.2 in 





Table 27. Beam Displacement Measurements for Shape Control Experiement 2 
(Nonlinear Model Voltages) 
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Actuator Inputs: 46.5/50/50/50 Volts 
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Figure 54. Beam Displacement Measurements for Shape Control Experiment 2 
(Nonlinear Model Voltages) (1,2 of 3) 
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Actuator Inputs: 141.1/150/150/150 Volts 
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Figure 55. Beam Displacement Measurements for Shape Control Experiment 2 
(Nonlinear Model Voltages) (3 of 3) 
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VI. CONCLUSIONS 


Piezoceramic actuators have been shown to provide an effective means of 
controlling the shape of a thin flexible beamlike structure. Experimental analysis has 
demonstrated that the Euler-Bernoulli model provides a good measure of the relative 
effectiveness of various actuator configurations and a reasonable model of actuator 
performance at low input voltage. However, the true nonlinear behavior of piezoceramic 
actuators and the effects of hysteresis and transverse stresses not accounted for by the 
Euler-Bernoulli model have been shown to have significant ramifications for the use of 
piezoceramics in shape control. Hysteresis resulted in residual beam deformation in excess 
of 20% of full deformation after the removal of the initial actuator input voltage for the 
Navy Type II piezoceramics tested, but following displacement measurements were found 
to be highly repeatable. Nonlinear voltage-strain behavior of the piezoceramic produced 
significantly higher displacements than predicted by the Euler-Bernoulli model at higher 
actuator input voltages. Transverse stresses were found to result in reduced effectiveness 
for a given thickness of actuator on wider beams due to the resulting increase in the 
moment of inertia of the beam cross-section. All of these effects result in variation of the 
relationship between actuator curvature and input voltage from the linear behavior 
predicted by the Euler-Bernoulli model. Precision shape control applications must 
therefore consider these effects and should employ a control system incorporating some 
form of strain and/or position feedback to ensure the desired deformation of the beam is 
achieved. 

The amount of beam deformation available from a piezoceramic actuator has been 
shown to depend on the dimensions and material properties of the structure and 
piezoceramic. High piezoceramic elastic modulus, piezoelectric coefficient and operating 
voltage are all desired actuator properties. The optimum total actuator thickness for a 
given beam/actuator configuration has been shown to be equal to a fixed percentage of the 
beam thickness, while the total beam curvature attainable for the optimum beam/actuator 


thickness ratio has been shown to be inversely related to the thickness of the beam. 
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The location of actuators on a beam structure is crucial to determining how closely 
a desired deformation profile can be approximated. Simultaneous optimization of the 
locations and input voltages of a fixed set of actuators to achieve a desired beam 
deformation profile proved to be unreliable due to the differences in the order of the 
actuator location and voltage terms in the optimization cost function. Using embedded 
Nelder and Mead simplex algorithms to separately optimize actuator locations and input 
voltages was found to produce much more reliable results, converging to the same 
optimum solution for a variety of initial conditions. Optimization of locations for different 
numbers of actuators making up a fixed total actuator length has shown that dividing a 
given actuator length into multiple small actuators allows a desired beam shape to be 
approximated more closely than can be achieved with a lesser number of longer actuators. 

Certain sets of initial conditions were found to cause the optimization algorithm to 
fail to converge or converge to an erroneous result. The algorithm was also found to be 
susceptible to converging to a local minimum of the cost function, when one existed and 
initial conditions near the local minimum were selected. Consideration must therefore be 
given to the presence of inflection points in the desired shape function, which can result in 
local minima of the cost function, when selecting initial actuator locations for the 
optimization algorithm. Knowledge of the approximate values of the optimum actuator 
locations should be used in selecting initial conditions whenever possible. Multiple runs of 
the optimization program should be performed using different sets of initial conditions to 
ensure that the globally optimum actuator locations are attained. 

Recommendations for further research include improving the mathematical model 
of the beam to consider transverse effects, incorporating constraints and techniques to 
avoid local minima in the optimization algorithm and exploring the use of feedback to 
control the shape of the beam. A plate theory model can be used to consider the effects of 
the curvature about the beam’s longitudinal axis caused by the actuator acting in both the 
longitudinal and transverse directions for beam’s with high width-to-thickness ratios. 
Modifying the optimization algorithm to include constraints would allow for the possibility 


of desired shape functions which can be most closely approximated by placing actuators 


126 


end to end or by actuators providing the maximum possible curvature. Techniques such 
as genetic algorithms or simulated annealing could also be included in the optimization 
algorithm to improve the probability of attaining the global minimum rather than a local 
minimum of the cost function. Multiple strain and/or displacement sensors could be 
utilized to provide feedback of the beam's shape to control inputs to the piezoceramic 
actuators. Sensors such as the laser displacement sensor used in the experiments 
conducted for this thesis and state-of-the-art fiber optic strain sensors are strong 


candidates for this application. 
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% PROGRAM TO COMPARE PERFORMANCE OF PIEZOCERAMIC ACTUATOR 
% CONFIGURATIONS 

% 

% This program calculates and plots the following data: 

% 1) Beam curvature per unit voltage 

% 2) Beam curvature at 20 Volts/mil of piezo thickness 

Vo 3) Beam surface strain per unit voltage (Side of beam 

% X Opposite actuator for unsymmetric configurations) 

% 4) Neutral axis distance from center of beam 

% versus actuator thickness or actuator width for piezoceramic 
y/o actuators bonded to a cantilever beam, given input properties 
% of the beam, actuators and the bond between the actuator and 
% the beam. 

Yo 

% The following actuator configurations are evaluated: 

% l)Asingle actuator bonded to one side of the beam 

% 2) A single actuator bonded to each side of the beam 

Vo 3) A dual (stacked) actuator bonded to one side of the beam 
% 4) A dual (stacked) actuator bonded to each side of the beam 


% Set number of points for plots 
Numpts = 100; 


% Set properties of beam 
% Material name 
Beammatl = '7075-T6 Aluminum’; 
% Width (m) 
\ = 0.0381; % 1.5 inches 
Vo Thickness 
t = 0.001575; 
% Elastic modulus (N/m‘%2) 
Е = 7.2е10; 


% Set properties of PZT actuator 

“% Material name 
Piezonam = “DOD Type II PZT’; 

% Constant width (m) for plots vs. actuator thickness 
wp = 0.0381; % 1.5 inches 

% Width range (m) for plots vs. actuator width 
wprange = linspace(0.1,1.0,Numpts)*w; 

% Constant thickness (m) for plots vs. actuator width 
tp = 0.00026; 
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% Thickness range (m) for plots vs. actuator thickness 
tprange = linspace(0.0001,0.002, Numpts); 

% Elastic modulus (N/m‘2) 
Ep = 6.1e10; 

% Piezoelectric coefficient (m/V) 
131 = 1.71е-10; % DOD Type II PZT 


Yo Set properties of bond 
% Width (m) 
wb = wp; 
% Thickness (m) 
tb = 0; 
?/o Elastic modulus (N/m^2) 
Eb = 1.78e9; 


% Set electric field limitation of piezoceramic (volts/mil) 
Volt_mil = 20; % Volts/mil 


?/o Calculate electric field limitation of piezoceramic (volts/m) 
V_tlimit = Volt_mil/0.0000254; % Volts/m = Volts/mil/(0.0000254m/mil) 


y/o Allow plots to be generated vs actuator width or thickness 
Plottype — input('Generate plots vs actuator thickness (t) or width (w)? ','s'); 


Yo Allow actuator thickness to be plotted as % of beam thickness 
% or actual thickness for plots vs actuator thickness 


if Plottype ~= 'w' 

Axistype = input('Plot actuator thickness as % of beam thickness (p) or actual thickness 
(а)? ',5); 
end 


Yo Create actuator width or thickness range in accordance 
% with plot type 
if Plottype = 'w' 
w = w*ones(1,Numpts); 
wp = wprange; 
wb = wp; 
else 
— t*ones(1,Numpts); 
tb = tb*ones(1,Numpts); 
tp = tprange; 
end 


% Calculate parameters for Crawley/Anderson formulation 
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Yo Single Actuator 

EA] = E*w.*t + Eb*wb.*tb + Ep*wp.*tp; 

ES1 = Eb*wb. *tb.*(t/2+tb/2) + Ep*wp.*tp.*(t/2+tb+tp/2); 

ЕЕ НАЗЛО Eb*wb. *tb *((t/2)^2 «-(t/2).*tb-- tb.^2/3)  ... 
Ep*wp* tp. *((t/2+tb).%2 + (t/2+tb).*tp + tp.%2/3); 

Р1 = Ep*d31*wp; 

MI - Ep*d31*wp.*(t/2 + tb + tp/2); 

EpsO 1 - (P1.*EII-MI.*ES1)./(EAI.*EII-ES1.^2), 

K V1 - MI /EII - ESI./EI1.*(P1.*EI1 - MI.*ES1).(EAI.*EII - ESI1.^2); 

yO 1=Eps0 1/K V1; 


% Symmetric Pair of Actuators 

EA2 = E*w.*t + 2*Eb*wb.*tb + 2*Ep*wp.*tp; 

DIPEE*w:*t283/12 --2*BEb*wbsttb.*((1/2).^2i8 (2). Mb tb.^2/8yere 
2*Ep* wp. *tp.*((t/2+tb).%2 + (t/2+tb).*tp + tp.%2/3); 

M2 = 2*Ep*d31*wp.*(t/2 + tb + tp/2); 

К У2 = М2./Е12; 


% Stacked Pair of Actuators 

EA12 = E*w.*t + 2*Eb*wb.*tb + 2*Ep*wp.*tp; 

ES12 = Eb*wb.*tb.*(t+2*tb + tp) + Ep*wp.*tp.*(t+3 *tb+2*tp); 

El12 = E*w.*t.43/12 + Eb*wb.*tb.*((t/2).%2 + (t/2+tbttp).%2 + ... 
(t+5/3*tb+tp).*tb) + Ep*wp.*tp.*((t/2+tb).%2 + ... 
(t/2+2*tb+tp).%2 + (t+3 *tb+5/3*tp).*tp); 

P12 = 2*Ep*d3 1*wp; 

M12 = Ep*d31*wp.*(t/2 + tb + tp/2 + t/2 + 2*tb + 3/2*tp); 

EpsO_12 = (P12.*EI12-M12.*ES12)./(EA12.*EI12-ES12.%2); 

K V12-7 MI2/EII2 - ES12/EI12.*(P12.*EI12 - MI2.*ES12)/ ... 

(EA12.*EI12 - ES12.^2); 
у0 12 = Ерѕ0 12./К У12; 


% Symmetric Pair of Stacked Actuators 

EA22 = E*w.*t + 4*Eb*wb.*tb + 4*Ep*wp.*tp; 

EI22 7 E*w.*t.^3/12 * 2*Eb*wb.*tb.*((t/2).^2 + (t/2+tbttp).%2 +... 
(t+5/3*tb+tp).*tb) + 2*Ep*wp.*tp.*((t/2+tb).2 +. ... 
(t/242*tb--tp).^2 + (t+3*tb+5/3*tp).*tp); 

M22 - 2*Ep*d31*wp.*(t/2 + tb + tp/2 + t/2 + 2*tb + 3/2*tp); 

К У22 = М22./Е122; 


% Calculate strain per unit actuator input voltage at beam surface 
% for each actuator configuration 

StrainV] - Eps0 1-K V1.*t/2; % Single Actuator 

StrainV2 = K_V2.*t/2; % Symmetric Pair of Actuators 
ЅігаіпУ12 = Ерѕ0 12+К V12.*t/2; % Stacked Pair of Actuators 
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StramV22=K V22.*t/2; % Symmetric Pair of Stacked Actuators 


% Plot results 
figure 
if Plottype — 'w' 
plot(wp./w*100,K Vl,'b--,wp./w*100,K V12,'m-.,wp./w*100,K V2,r-, ... 
wp./w*100,K V22,'g.") 
text(0.15,0.6,[num2str(tp* 1000) ' mm Thick Actuator/ num2str(tb(1)*1000) ... 
' mm Bond Thickness |,'Units',"Normalized") 
xlabel(‘Actuator Width (%)') 
elseif Axistype — 'p' 
plot(tp./t*100,K V1,'b--'tp./t*100,K V12,'m-.,tp./t*100,K V2,r-' ... 
tp./t*100,K V22,g.) 
text(0.15,0.6,[num2str(wp/w*100) '9$? Width Actuator/' num2str(tb(1)* 1000) ... 
' mm Bond Thickness'],'Units','Normalized") 
xlabel('Actuator Thickness (% of Beam Thickness)') 
else 
plot(tp*1000,K V1,b--,tp*1000,K V12,m-.,tp*1000,K V2,r-, ... 
tp*1000,K V22,g.) 
text(0.15,0.6,[num2str(wp/w* 100) '% Width Actuator/ num2str(tb(1)*1000) ... 
' mm Bond Thickness ],'Units',"'Normalized') 
xlabel('Actuator Thickness (mm)') 
end 
ylabel('(rad/m/V)) 
title(‘Beam Curvature per Unit Voltage’) 
legend('Single','l Stacked','Sym. Pair','Sym. Stacked',- 1) 
text(0.15,0.9, Piezonam,'Units','Normalized') 
text(0.15,0.8,'Ep ^" num2str(Ep) ' N/m^2 431 = ' num2str(d31) ' m/ V'], ... 
"Units',"Normalized') 
text(0.15,0.7,[num2str(t(1)*1000) ' mm Thick ' Beammatl ' Beam'],'Units,, ... 
'Normalized') 
grid 


figure 
if Plottype — 'w' 
plot(wp./w*100,V tlimit*tp.*K V1,'b--,wp./w*100,V tlimit*tp.*K V12,'m-., ... 
wp./w*100,V tlimit*tp.*K V2,'r-wp./w*100,V tlimit*tp.*K V22,'g. 
text(0.15,0.6,[num2str(tp* 1000) ' mm Thick Actuator/ num2str(tb(1)*1000) ... 
' mm Bond Thickness'],'Units',"Normalized") 
xlabel(‘Actuator Width (“%)') 
elseif Axistype = 'p' 
plot(tp.*100,V_tlimit*tp.*K_V1,'b--',tp./t*100,V_thimit*tp.*K_V12,'m-"', ... 
tp./t*100,V_tlimit*tp.*K_V2,'r-',tp./t*100,V_tlimit*tp.*K_V22,'g.') 
text(0.15,0.6,[num2str(wp/w* 100) '% Width Actuator/ num2str(tb(1)* 1000) ... 
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: mm Bond Thickness ],'Units',"Normalized") 
xlabel('Actuator Thickness (Vo of Beam Thickness) 
else 
plot(tp*1000,V_tlimit*tp.*K_V1,'b--',tp* 1000, V_tlimit*tp.*K_ V12,'m-.", ... 
tp*1000,V tlimit*tp.*K V2,r-,tp*1000,V tlimit*tp.*K V22,g.') 
text(0.15,0.6,[num2str(wp/w* 100) '% Width Actuator/' num2str(tb(1)*1000) ... 
'mm Bond Thickness'],"Units',"Normalized') 
xlabel('Actuator Thickness (mm) 
end 
ylabel('‘Curvature (rad/m)') 
title(['Beam Curvature at ' num2str(Volt mil) ' Volts/mil']) 
legend('Single','1 Stacked','Sym. Pair','Sym. Stacked',-1) 
text(0.15,0.9, Piezonam,'Units','Normalized') 
text(0.15,0.8,['Ep — ' num2str(Ep) ' N/m^2 431 -' num2str(d31) ' m/V'|, ... 
'Units', Normalized’) 
text(0.15,0.7,[num2str(t(1)* 1000) ' mm Thick ' Beammatl ' Beam'],'Units,, ... 
'Normalized") 
grid 


figure 
if Plottype == 'w' 
plot(wp./w* 100,StrainV1* 1e6,'b--',wp./w* 100,StrainV 12*1e6,'m-.’, ... 
wp./w*100,StrainV2* 1e6,'r-',wp./w* 100,StrainV22* 1e6,'g.') 
text(0.15,0.6,[num2str(tp* 1000) ' mm Thick Actuator/ num2str(tb(1)*1000) ... 
' mm Bond Thickness ],'Units','Normalized') 
xlabel('Actuator Width (96)') 
elseif Axistype == 'р' 
plot(tp./t* 100, StrainV 1 *1e6,'b--' tp./t* 100, StrainV 12*1e6,'m-.", ... 
tp./t* 100, StrainV2* 1e6,'r-tp./t* 100, StrainV22* 1e6,'g.") 


text(0.15,0.6,[num2str(wp/w* 100) '% Width Actuator/ num2str(tb(1)*1000) ... 


' mm Bond Thickness ],'Units',"Normalized') 
xlabel('Actuator Thickness (9o of Beam Thickness) 
else 
plot(tp* 1000,StrainV 1*1e6,'b--' tp* 1000,StrainV 12*1e6,'m-.', ... 
tp* 1000, StrainV2*1e6,'r-' tp* 1000, StrainV22* 1e6,'g.") 


text(0.15,0.6, [num2str(wp/w* 100) '9s Width Actuator/ num2str(tb(1)*1000) ... 


‘mm Bond Thickness ],'Units',' Normalized") 

xlabel('Actuator Thickness (mm)') 
end 
ylabel('Strain/Voltage (microstrain/V)') 
title" Surface Axial Strain per Unit Actuator Voltage") 
legend('Single','1 Stacked','Sym. Pair','Sym. Stacked',-1) 
text(0.15,0.9 Piezonam, Units',"'Normalized") 
text(0.15,0.8,['Ep 7 ' num2str(Ep) ' N/m^2  d31 -' num2str(d31) ' m/V', ... 
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'Units';'Normalized") 

text(0.15,0.7,[num2str(t(1)* 1000) ' mm Thick ' Beammatl ' Beam'],'Units,, ... 
'Normalized') 

grid 


figure 
if Plottype — 'w' 
plot(wp./w*100,y0 1*1000,;b--'wp./w*100,y0 12*1000,'m-.") 
text(0.15,0.6,[num2str(tp* 1000) ' mm Thick Actuator/' num2str(tb(1)*1000) ... 
' mm Bond Thickness'|,'Units',"'Normalized") 
xlabel(‘Actuator Width (%)') 
ylabel('yO (mm) 
elseif Axistype — 'p' 
plot(tp./t*100,y0 1./t*100,'b--'tp./t*100,y0 12/t*100,'m-.') 
text(0.15,0.6,[num2str(wp/w* 100) '9; Width Actuator/' num2str(tb(1)*1000) ... 
' mm Bond Thickness'],"Units',"Normalized') 
xlabel('Actuator Thickness (vo of Beam Thickness)') 
ylabel('yO (Yo of Beam Thickness)') 
else 
plot(tp*1000,y0 1*1000,b--'tp*1000,y0 12*1000,'m-.") 


text(0.15,0.6,[num2str(wp/w* 100) '% Width Actuator/’ num2str(tb(1)* 1000) ... 


' mm Bond Thickness'],'Units',"Normalized') 
xlabel('Actuator Thickness (mm)') 
ylabel(‘yO (mm)') 
end 
title("Neutral Axis Distance from Center of Beam’) 
legend('Single','1 Stacked',-1) 
text(0.15,0.9, Piezonam, Units','Normalized') 
text(0.15,0.8,['Ep ^ ' num2str(Ep) ' N/m^2 431 = ' num2str(d31) ' m/V!], ... 
'Units','Normalized") 
text(0.15,0.7,[num2str(t(1)*1000) ' mm Thick ' Beammatl ' Beam'],'Units’, ... 
'Normalized') 
grid 
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% PROGRAM TO DETERMINE OPTIMUM ACTUATOR LOCATIONS 
% 

Yo This program determines the optimum actuator locations and input 

% voltages to achieve the minimum error between the actuator deformed 

% beam shape and a desired beam shape function for a cantilever 

y/o beam with piezoceramic actuators of specified length, width, thickness 
Yo and material properties bonded to both sides of the beam. 


% Set desired shape function and coefficients 
Shapfunc = 1 % у = C*x^2 
% Shapfunc = 2 % y = C*(1-cos(mx)) 


C = 0.001; 
% C = 0.0005; 
m = pi, 


y/o Set beam properties 


Г = 1, % Веат length (m) 

уу = 0.0508; % Beam width (m) 

% t = 0.001575;  . Yo Beam thickness (m) 

t — 0.001575/2 ?/o Beam thickness (m) 

E = 7.2e10; % Beam modulus of elasticity (N/m^2) 
% Set piezoceramic actuator properties 

% Navy Type Ii PZT 

n=4 % Number of actuators 

la = 0.1*ones(n,1); % Actuator lengths (m) 

wp = 0.0381; % Piezoceramic actuator width (m) 

tp 0.00027; % Piezoceramic actuator thickness (m) 
Ep = 6.1e10; % Piezoceramic actuator modulus of elasticity (N/m^2) 
431 = 1.71е-10; % Piezoelectric constant (m/V) 

Vmax = 150; % Maximum actuator input voltage (V) 
Vmin = -150; % Minimum actuator input voltage (V) 


% Enter properties of bond between piezoceramic actuator and beam 


wb = wp; % Bond width (m) (equal to actuator width) 

tb = 0; % Bond thickness (m) 

Eb = 1.78e9; % Bond modulus of elasticity (N/m^2) 

% Calculate properties of beam cross-section 

area — w.*t; % Beam cross-sectional area (m^2) 

areab — wb*tb; Yo Bond cross-sectional area (m^2) 

areap — wp*tp; % Piezoceramic actuator cross-sectional area (m^2) 


по! = w.*t.^3/12; 9 Beam cross-section moment of inertia (m^4) 
% Bond cross-section moment of inertia (m^4) 
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moib — wb*tb*(t^2/2 + t*tb + 2*tb%2/3); 

% Piezoceramic actuator cross-section moment of inertia (m4) 

moip = wp*tp*((t+2*tb)%2/2 + (t+2*tb)*tp + 2*tp^2/3), 

EI = E*moi+Eb*moib+Ep*moip; % Beam stiffness 

К = Ep*d31*wp*(t+2*tb+tp)/EI; % (Moment per unit voltage)/(Beam stiffness) 
K = K*ones(n, 1); 


% Set options for unconstrained voltage optimization algorithm 
Options = foptions’; 
Options(14) = 2000; % Maximum number of iterations (0 -> 100*n) 


% Set initial actuator locations for optimization algorithm 
x0 = [0.05 0.25 0.5 0.75] % 4 Actuators 


% Determine optimum voltage for initial actuator positions 
if Shapfunc == 
VO = fmins(‘beamerVu',zeros(n, 1), Options,[],x0,la,K,C,L,n) 
% VO = zeros(n,1) 
else 
V0 = fmins(‘bercos Vu',zeros(n, 1 ),Options,|[ ],x0,la,K,C,m,L,n) 
end 


% Calculate error cost function for initial actuator positions 
if Shapfunc = 1 
JO = beamerVu(V0,x0,la,K,C,L,n) 
else 
JO = bercosVu(V0,x0,la,K,C,m,L,n) 
end 


jo Determine optimum actuator locations 
if Shapfunc == 1 
flops(0) % Reset floating point operations counter 
Optx = fmins(‘beamerxu',x0, Options, [], V0,la,K,C,L,n, Options) 
Ops = flops % Output number of floating point operations 
else 
flops(0) % Reset floating point operations counter 
% Unconstrained locations 
Optx = fmins(‘bercosxu',x0,Options,[], V0,la,K,C,m,L,n,Options) 
Ops = flops % Output number of floating point operations 
end 


% Determine optimum input voltages for optimum actuator locations 
if Shapfunc = | 
OptV = fmins(beamerVu', VO Options,[], Optx,la, K, C,L,n) 
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else 
Opt V = fmins(‘bercosVu', V0,Options,[],Optx,la,K,C,m,L,n) 
end 


% Calculate error cost function for optimum actuator locations 
if Shapfunc == 1 
J = beamerVu(OptV,Optx.la,K,C,L,n) 
else 
J = bercosVu(OptV,Optx,la,K,C,m,L,n) 
end 


% Plot results 
% Set plot x range along beam length 
x = linspace(0,L, 100); 
% Calculate beam displacement for optimum actuator locations and voltages 
y = beamdisp(x,Optx,OptV,la,n,K); 
% Calculate beam displacement at actuator actuator endpoints 
yact = beamdisp({Optx;Optx+tla],Optx,OptV,la,n,K); 
% Calculate desired beam displacement along beam length 
if Shapfunc — 1 
yf = C*x.^2; 
else 
yf = C*(ones(size(x))-cos(m*x)); 
end 
figure 
plot(x,y,'r-',x, yf,'g.',[Optx;Optx+la], yact,'r*') 
legend(‘Actual','Desired’,-1) 
grid 
title('Actual vs. Desired Beam Shapes’) 
xlabel('X (m)') 
ylabel(‘Beam Displacement Y (m)') 
ifn == ] 
text(0.1,0.9,['Optimum Values for 1 Actuator'],"Units','Normalized") 
else 
text(0.1,0.9,['Optimum Values for ' int2str(n) ' Actuators'],'Units,, ... 








'Normalized") 
end 
ifn—1 
text(0.1,0.8;['x] = ' num2str(Optx(1)) ' m / V1 -' num2str(OptV(1)) ... 
' V','Units,'Normalized') 
text(0.1,0.7,['Error = ' num2str(J)],'Units','Normalized") 
elseif n — 2 
text(0.1,0.8,['x] =' num2str(Optx(1))'m / V1 =' num2str(OptV(1)) ... 
''V'],'Units',"Normalized’) 
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text(0.1,0.7,['x2 ^ ' num2str(Optx(2)) ' m / V2 -' num2str(OptV(2)) ... 
' V'],'Units','Normalized") 
text(0.1,0.6,['Error = ' nun2str(J)],'Units','Normalized') 
elseif n — 3 
text(0.1,0.8,['x1 —' num2str(Optx(1)) ' m / V1 ^' num2str(OptV(1)) ... 





' V'|,'Units','Normalized") 

text(0.1,0.7,['x2 ^ ' num2str(Optx(2)) ' m / У2 = ' num2str(OptV(2)) ... 
' V'],'Units','Normalized") 

text(0.1,0.6,['x3 — ' num2str(Optx(3)) ' m / V3 -' num2str(OptV(3)) ... 
' V'] 'Untts', Normalized') 


text(0.1,0.5,['Error = ' num2str(J)],'Umts' Normalized') 
elseif n — 4 

text(0.1,0.8,['x1 — ' num2str(Optx(1)) ' m / V1 -' num2str(OptV(1)) ... 
' V'],'Units',"Normalized') 

text(0.1,0.7,['x2 =' num2str(Optx(2))'m / V2 =' num2str(OptV(2)) ... 
' V''Units','Normalized') 

text(0.1,0.6,['x3 = ' num2str(Optx(3)) ' m / V3 7 ' num2str(OptV(3)) ... 
' V','Units','Normalized) 

text(0.1,0.5,['x4 = ' num2str(Optx(4)) ' m / У4 = ' num2str(OptV(4)) ... 
' V','Units','Normalized") 

text(0. 1,0.4,['Error = ' num2str(J)],"Units',"Normalized’) 

else % (п = 5) 

{ех((0.1,0.8, ['х1 = ' num2str(Optx(1)) ' m / УІ = ' num2str(OptV(1)) ... 
' V'|,'Units',"Normalized') 

text(0.1,0.7,['x2 7 ' num2str(Optx(2)) ' m / V2 7 ' num2str(OptV(2)) ... 
' V'|,'Units','Normalized") 

text(0.1,0.6,['x3 ^ ' num2str(Optx(3)) ' m / V3 -' num2str(OptV(3)) ... 
' V'|,'Units',"Normalized") 

text(0.1,0.5,['x4 7 ' num2str(Optx(4)) ' m / V4 - ' num2str(OptV(4)) ... 
' V'],'Units','Normalized") 


text(0.1,0.4,['x5 ^ ' num2str(Optx(5)) ' m / V5 =' num2str(OptV(5)) ... 


' V],'Units','Normalized^) 
text(0.1,0.3,['Error 7 ' num2str(J) |, Units',"Normalized') 
end 


Yo Allow data to be saved 

Savedata — input('Save data? (y/n) ','s); 

if Savedata — 'y' 
Savefile = input('Enter filename (without .MAT extension) ','s'); 
eval(['save ' Savefile]) 

end 
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% FUNCTION TO DETERMINE ERROR BETWEEN ACTUAL BEAM SHAPE AND 
% DESIRED PARABOLIC SHAPE FOR UNCONSTRAINED ACTUATOR 
% LOCATION OPTIMIZATION 

% 

% function f = beamerxu(x, V0,la,K,C,L,n,Options) 

% 

Yo This function calculates the value of the cost function 'f 

% for a cantilever beam of length 'L' with 'n' attached 

% piezoceramic actuators at locations 'x' with lengths 'la' and 

?/o curvature per unit input voltage 'K' and a desired parabolic 

Yo beam shape function y ^ C*x^2. 


function f 2 beamerxu(x, VO,la, K C,L,n,Options) 


% Determine optimum actuator input voltages for input actuator 
% locations 'x' 
V = fmins(‘beamer Vu', VO, Options, | ],x,la,K,C,L,n); 


% Calculate summations for use in calculating beam shape 

% polynomial coefficients 

sumK VI = zeros(n, 1); 

sumK VIxl = zeros(n, 1); 

torr — in 
sumK Vl(i) » sum(K(1:1).*V(1:).*la(1:))); % Sum of Ki*Vi*li 
sumKVlxl(1) » sum(K(1:1).*V(1:1).*1a(1:1). *(x(1:1)*1a(1:1)/2)); 

99 Sum of Ki*Vi*li*(xi*l/2) 
end 


% Calculate beam shape polynomial coefficients 
a2 = K.*V/2; 

al = -K.*V.*x; 

al(2:n) = al(2:n) + sumK VI1(1:n-1); 

a0 7 K.*V.*x.^2/2; 

a0(2:n) = a0(2:n) - sumKVlIxl(1:n-1); 

bl = sumK V1; 

bO = -sumK VIxl; 


% Create nx1 vector containing value of 'C' 
Cvect = C*ones(n, 1); 


% Calculate value of cost function 

i: C^2*x(1)^5/5: 

f= f+ sum((a2-Cvect).%2.*((x+la).%5 - x.^5)y/5; 
f= f+ sum((a2-Cvect).*al .*((xt+la).%4 - x.%4))/2; 
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f = f + sum((2*(a2-Cvect).*a0+al .%2).*((x+la).%3 - x.^3))/3; 

f= f+ sum(al .*a0.*((xtla).%2 - x.%2)); 

f ^ f  sum(a0.^2.*la); 

f - f* C^2*sum(x(2:n).^5 - (x(1:n-1)-Ha(1:n-1)).^5y/5; 

f= f - C*sum(bl(1:n-1).*(x(2:n).^4 - (x(1:n-1)*a(1:n-1)).^4))/2; 

f= f+ sum((-2*C*b0(1:n-1)*b1(1:n-1).^2).*(x(2:n).^3 - ... 
(x(1:n-1)+la(1:n-1)).%3))/3; 

f—f-* sum(bl(1:n-1).*bO(1:n-1).*(x(2:n).^2 - (x(1:n-1)-a(1:n-1)).^2)); 

f= f+ sum(b0(1:n-1).%2.*(x(2:n) - x(1:n-1)-la(1:n-1))); 

f= f + C*2*(L%5-(x(n)+la(n))*5)/5; 

f= f - C*b1(n)*(L“4-(x(n)+la(n))4)/2; 

f—-f- (-2*C*bO(n)*bl(n)^2)*(L^3-(x(n)-a(n))^3)/3; 

f ^ f * bl(n)*bO(n)*(L^2-(x(n)-la(n))^2); 

f= f+ bO(n)*2*(L-x(n)-la(n)); 


% FUNCTION TO DETERMINE ERROR BETWEEN ACTUAL BEAM SHAPE AND 
% DESIRED PARABOLIC SHAPE FOR UNCONSTRAINED ACTUATOR 
% VOLTAGE OPTIMIZATION 

Yo 

% function f = BeamerVu(V,x,la,K, C,L,n) 

% 

% This function calculates the values of the cost function 'Р 

Yo for a cantilever beam of 

% length 'L' with 'n' attached piezoceramic actuators and a 

% desired parabolic beam shape function y = C*x”%2, where the 

% elements of the input vector 'V' are the actuator input 

% voltages, the elements of 'x' are the distances 

% of the actuators from the root of the beam, and the 

y/o elements of 'la' are the actuator lengths. 


function f= BeamerVu(V,x,la,K,C,L,n) 


% x = Actuator locations (m) 
% V = Actuator voltages (V) 
% la = Actuator lengths (m) 


% Calculate summations for use in calculating beam shape 

% polynomial coefficients 

sumK V1 = zeros(n, 1 ); 

sumK V1xl = zeros(n, 1); 

for i= l:n 
sumKVl(i) ? sum(K(1:1).*V(1:1).*1a(1:1)); 96 Sum of Ki* Vi*li 
sumKVlxl(1) ^? sum(K(1:1). *V(1:1).*1a(1:1).* (x(1:1)Ha(1:1)/2)); 
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% Sum of Ki* Vi*li*(xitli/2) 
end 


% Calculate beam shape polynomial coefficients 
a2 = K.*V/2; 

al = -K.*V.*x; 

al(2:n) = а1(2:п) + sumKVI(1:n-1); 

a0 — K.*V.*x.^2/2; 

a0(2:n) 7 a0(2:n) - sumKVIxl(1:n-1); 

bl = sumK V1; 

bO = -sumK VIxl; 


% Create nx1 vector containing value of 'C' 
Cvect — C*ones(n,1); 


% Calculate value of cost function 

ИЕ 672 +х(1)2575; 

f= f + sum((a2-Cvect).%2.*((xtla).%5 - x.%5))/5; 

f= f+ sum((a2-Cvect).*al.*((xtla).%4 - x.^4))/2; 

f= f + sum((2*(a2-Cvect).*a0*a1.^2).*((x-Ha).^3 - x.^3))/3; 

f= f+ sum(al.*a0.*((x+la).%2 - x.^2)); 

f= f+ sum(a0.%2.*la); 

f= f+ C^2*sum(x(2:n).^5 - (x(1:n-1)-a(1:n-1)).^5)/5; 

f — f- C*sum(bl1(1:n-1).*(x(2:n).^4 - (x(1:n-1)-a(1:n- n. А 

f—f-* sum((-2*C*bO(1:n-1)*b1(1:n-1).^2).*(x(2:n).^3 - ... 
(x(1:n-1)+la(1:n-1)).%3))/3; 

f— f * sum(bl1(1:n-1).*bO(1:n-1).*(x(2:n).^2 - (x(1:n-1)-Ha(1:n-1)).^2)); 

Ё = f * sum(b0(1:n-1).^2.*(x(2:n) - x(1:n-1)-la(1:n-1))); 

f—-f-*C^2*(L^5-(x(n)-a(n))^5)/5; 

f —-f- C*bl(n)*(L^4-(x(n)-la(n))^4)2; 

f—-f-(-2*C*bO(n)bl(ny)2)*(L^3-(x(n)*la(n))^3)y3; 

f= f + bl(n)*bO(n)*(L^2-(x(n)-a(n))^2); 

f= f+ b0(n)^2*(L-x(n)-la(n)); 
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Yo FUNCTION TO DETERMINE ERROR BETWEEN ACTUAL BEAM SHAPE AND 
% DESIRED C*(1-cos(mx)) SHAPE FOR UNCONSTRAINED ACTUATOR 
% LOCATION OPTIMIZATION 

% 

% function f = bercosxu(x, VO,la,K,C,m,L,n, Options) 

% 

Yo This function calculates the value of the cost function 'f 

% for a cantilever beam of length 'L' with 'n' attached 

% piezoceramic actuators at locations 'x' with lengths 'la' and 

% curvature per unit input voltage 'K' and a desired sinusoidal 

% beam shape function y = C*(1-cos(mx)). 


function f = bercosxu(x, V0,la,K,C,m,L,n, Options) 


% Determine optimum actuator input voltages for input actuator 
y/o locations 'x' 
V — fmins('bercosVu', VO Options,[],x,la, K, C,m,L,n); 


% Calculate summations for use in calculating beam shape 

% polynomial coefficients 

sumK VI = zeros(n, 1); 

sumK VIxl = zeros(n, 1); 

fori- Ln 
sumKVl(1) = sum(K(1:1).*V(1:1).*1a(1:1)); 9o Sum of Ki*Vi*li 
sumKVlxl(i) » sum(K(1:1). *V(1:1).*la(1:1). *(x(1:1)*a( 1:1)/2)); 

9o Sum of Ki*Vi*li*(xi-li/2) 
end 


% Calculate beam shape polynomial coefficients 
a2 — K.*V/2; 

a] = -K.*V.*x, 

al(2:n) = al(2:n) + sumK VI(1:n-1); 

a0 — K.*V.*x.^2/2; 

a0(2:n) 7 a0(2:n) - sumK Vlxl(1:n-1); 

b1 = sumK VI; 

b0 = -sumK Vlxl; 


% Create nx1 vector containing value of 'C' 
Cvect = C*ones(n,1); 


% Calculate value of cost function 

f= C’2*(x(1)-2/m* sin(m*x(1)) + L/2 + 1/4/m*sin(2*m*L)); 
f=f+ sum(a2.%2.*((xtla).%5 - x.%5))/5; 

f= f+ sum(a2.*al.*((xtla).%4 - x.%4))/2; 
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f = f * sum((2*a2.*(a0-Cvect)*al.^2).*((x-a).^3 - x.^3))/3; 

f —f- sum(al.*(a0-Cvect).*((x-a).^2 - x.^2)); 

f= f+ sum((a0-Cvect).%2.*la); 

f = f + 2*C*sum(a2.*(2/m^2*((x-Ha).*cos(m*(x-Ha))-x.* cos(m*x)))); 

f = f + 2*C*sum(a2.*(((x+la).^2/m-2/m^3 *ones(n,1)).*sin(m*(x+a))- ... 
(x.^2/m-2/m^3 *ones(n, 1)).*sin(m*x))); 

f= f+ 2*C*sum(al.*((cos(m*(x+la))-cos(m*x))/m^2 +... 
((x+la).*sin(m*(x+a))-x.*sin(m*x))/m)); 

f = f + 2/m*C*sum((a0-Cvect). *(sin(m*(x+la))-sin(m*x))); 

f — f * sum(bl(1:n-1).^2.*(x(2:n).^3 - (ха п-1)На(1:п-1)).73))/3; 

Ё= {+ sum(bl(1:n-1).*(bO(1:n-1)-Cvect(1:n-1)).*(x(2:n).^2 - ... 
(x(1:n-1)-a(1:n-1)).^2)); 

f — f * sum((b0(1:n-1)-Cvect(1:n-1)).^2.*(x(2:n) - x(1:n-1)-la(1:n-1))); 

f—f-*2*C*sum(bl(1:n-1).*(cos(m*x(2:n))-cos(m*(x(1:n-1)-a(1:n-1)))))/m^2; 

f—f-* 2*C*sum(bl(1:n-1).*(x(2:n).*sin(m*x(2:n))-(x(1:n-1)-a(1:n-1)).* ... 
sin(m*(x(1:n-1)+la(1:n-1)))))/m; 

f= f+ 2*C/m*sum((b0(1:n-1)-Cvect(1:n-1)).*(sin(m*x(2:n))- ... 
sin(m*(x(1:n-1)-a(1:n-1))))); 

f—f-*bl(n).^2.*(L.^3 - (x(n)a(n)).^3)/3; 

f— f- bl(n).*(bO(n)-C).*(L.^2 - (x(n)-Ha(n)).^2); 

f—f- (bO0(n)-C).^2.*(L - x(n)-la(n)); 

f= f + 2*C*bl(n).*(cos(m*L)-cos(m*(x(n)-a(n)))/m^2; 

f—-f-2*C*bl(n).*(L.*sin(m*L)-(x(n)-la(n)).*sin(m*(x(n)-a(n))))/m; 

f = f+ 2*C/m*(b0(n)-C).*(sin(m*L)-sin(m*(x(n)+la(n)))); 


% FUNCTION TO DETERMINE ERROR BETWEEN ACTUAL BEAM SHAPE AND 
% DESIRED C*(1-cos(mx)) SHAPE FOR UNCONSTRAINED ACTUATOR 
% VOLTAGE OPTIMIZATION 

% 

% function f = bercosVu(x, VO,la,K,C,m,L,n, Options) 

% 

Vo This function calculates the value of the cost function '£ 

Yo for a cantilever beam of length 'L' with 'n' attached 

% piezoceramic actuators at locations 'x' with lengths 'la' and 

Yo curvature per unit input voltage 'K' and a desired sinusoidal 

% beam shape function y = C*(1-cos(mx)). 


function f= bercosVu(V,x,la,K,C,m,L,n,Options) 
% x = Actuator locations (m) 


% V = Actuator voltages (V) 
Yo la = Actuator lengths (m) 
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% Calculate summations for use in calculating beam shape 

% polynomial coefficients 

sumK VI = zeros(n,1); 

sumK VIxl = zeros(n, 1); 

fori- ln 
sumKVl(1) » sum(K(1:1).*V(1:1).*la(1:1)); 9 Sum of Ki*Vi*li 
sumKVlxl(1) 7 sum(K( 1:1). *V(1:1). *la(1:1).*(x(1:1)a(1:1)/2)); 

% Sum of Ki* Vi*li*(xiHi/2) 
end 


% Calculate beam shape polynomial coefficients 
a2 = K.*V/2; 

al = -K.*V.*x; 

а1(2:п) = а1(2:п) + sumKVI(1:n- 1), 

а0 = К.*У.*х.^2/2; 

a0(2:n) 7 a0(2:n) - sumKVIxl(1:n-1); 

bl = sumK VI; 

bO = -sumK VIx1; 


% Create nx] vector containing value of 'C' 
Cvect = C*ones(n, 1); 


% Calculate value of cost function 

f= C*2*(x(1)-2/m* sin(m*x(1)) + L/2 + 1/4/m*sin(2*m*L)); 

f= f+ sum(a2.%2.*((xta).%5 - x.^5)/5; 

f= f+ sum(a2.*al.*((xtla).%4 - x.%4))/2; 

f= f+ sum((2*a2.*(a0-Cvect)+a1.%2).*((xtla).%3 - x.%3))/3; 

f= f+ sum(al .*(a0-Cvect).*((xtla).%2 - x.^2)); 

f = f + sum((a0-Cvect).”2.*la); 

f= f+ 2*C*sum(a2.*(2/m^2*((x-Ha).* cos(m*(x-la))-x. *cos(m*x)))); 

f—-f-*2*C*sum(a2.*(((x*la).^2/m-2/m^3 *ones(n, 1)). *sin(m*(x-a))- ... 
(x.^2/m-2/m^3 * ones(n,1)).*sin(m*x))); 

f-f-*2*C*sum(al.*((cos(m*(x-*la))-cos(m*x))/nm72 +... 
((x+la).*sin(m*(x+la))-x.*sin(m*x))/m)); 

f= f+ 2/m*C*sum((a0-Cvect). *(sin(m*(x+la))-sin(m*x))); 

f= f+ sum(b1(1:n-1).%2.*(x(2:n).43 - (x(1:n-1)+la(1:n-1)).%3))/3; 

f= f+ sum(b1(1:n-1).*(b0(1:n-1)-Cvect(1:n-1)).*(x(2:n).%2 - ... 
(x(1:n-1)+la(1:n-1)).%2)); 

f= f + sum((b0(1:n-1)-Cvect(1 :n-1)).%2.*(x(2:n) - x(1:n-1)-la(1:n-1))); 

f= f+ 2*C*sum(b1(1:n-1).*(cos(m*x(2:n))-cos(m*(x(1:n-1)+a(1:n-1)))))/m“2; 

f= f+ 2*C*sum(b1(1:n-1).*(x(2:n).*sin(m*x(2:n))-(x(1 :n-1)+la(1:n-1)).* ... 
sin(m*(x(1:n-1)+la(1:n-1)))))/m; 

f = f + 2*C/m*sum((b0(1:n-1)-Cvect(1 :n-1)).*(sin(m*x(2:n))- ... 
sin(m*(x(1:n-1)+a(1:n-1))))); 
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f—-f-tbl(n)^2.*(L.^3 - (x(n)-a(n)).^3y3; 

f— f bl(n).*(bO(n)-C).*(L.^2 - (x(n)a(n)).^2), 

f— f t (b0(n)-C).^2.*(L - x(n)-la(n)); 

f= f+ 2*C*b1(n).*(cos(m*L)-cos(m*(x(n)+la(n))))/m”2; 

f= f+ 2*C*bl1(n).*(L.*sin(m*L)-(x(n)+la(n)). *sin(m*(x(n)+la(n))))/m; 
f= f+ 2*C/m*(b0(n)-C).*(sin(m*L)-sin(m*(x(n)+la(n)))); 


% FUNCTION TO DETERMINE THE DEFLECTED SHAPE OF A CANTILEVER 
% BEAM DUE TO ATTACHED PIEZOCERAMIC ACTUATORS 
% 

Yo function y = beamdisp(x,xa,v,la,n,K) 

Yo 

Yo This function returns a vector of transverse displacements, y, 

% of a cantilever beam with 'n' attached piezoceramic actuators, 

% given the starting locations of the actuators, 'xa', the lengths 

% of the actuators, 'la', the actuator input voltages, 'v', and the 

% vector 'K' containing the actuator curvatures per unit voltage. 

Yo 'X 1s a Vector containing a range of axial distances from the 

y/o root of the beam at which 'y' is to be evaluated and 'xa', у 

% апа 'la' are all nx1 vectors. 


function y = beamdisp(x,xa,v,la,n,K) 
y = zeros(size(x)); 
for m = 1:length(x) 
if x(m) «— xa(1) 6x prior to start of actuator 
y(m) = 0; 
elseif xa(1) < x(m) & x(m) «- xa(1)-a(1) 
y(m) - K(1)*v(1)2*(x(m)-xa(1))^2; 
elseif x(m) <= xa(n)+la(n) 
for k= 1:n 
if xa(k) <= x(m) & x(m) <= xa(k)+la(k) 
y(m) » K(k)*v(ky2*(x(m)-xa(k))^2*sum(K(1:k-1).*v(1:k-1).* ... 
la(1:k-1).*(x(m)*ones(size(xa( 1:k-1)))-xa( 1:k-1)-la(1:k-1)/2)); 
elseif k<n 
if xa(k)+la(k) <= x(m) & x(m) <= xa(k+1) 
y(m) = sum(K(1:k).*v(1:k).*la(1:k).*(x(m)* ones(size(xa(1:k)))- ... 
xa(1:k)-la(1:k)/2)); 
end 
end 
end 
else 
y(m) 7 sum(K(1:n).*v(1:n).*la(1:n).*(x(m)*ones(size(xa(1:n)))- ... 
xa(1:n)-la(1:ny2)), 
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end 
end 


% PROGRAM TO DETERMINE OPTIMUM ACTUATOR LOCATIONS 
% FROM A GIVEN DISCRETE SET OF POSSIBLE LOCATIONS 

% 

% This program determines the optimum actuator location and input 

% voltage from a specified discrete range of possible actuator 

% locations to achieve the minimum error between the actuator deformed 
% beam shape and a desired parabolic shape function for a cantilever 

% beam with piezoceramic actuators bonded to both sides of the beam. 

% 1, 2 or 3 actuators can be specified 


% Coefficient of desired parabolic shape function 


С = 0.001; 

% Enter beam properties 

L=1; Yo Beam length (m) 

w — 0.0508; % Beam width (m) 

t = 0.001575; 7o Beam thickness (m) 

E = 7.2e10; % Beam modulus of elasticity (N/m“2) 
% Enter piezoceramic actuator properties 

% Navy Type H PZT 

n= 3; % Number of actuators 

la = 0.1*ones(n,1); % Actuator lengths (m) 

wp = 0.0381; % Piezoceramic actuator width (m) 

tp = 0.00027; % Piezoceramic actuator thickness (m) 
Ep = 6.le10; % Piezoceramic actuator modulus of elasticity (N/m^2) 


131 = 1.71е-10; 96 Piezoelectric constant (m/V) 


o Enter properties of bond between piezoceramic actuator and beam 
% Enter propert f bond bet p tuat db 


wb = wp; % Bond width (m) (equal to actuator width) 

tb = 0; % Bond thickness (m) 

Eb = 1.78e9; % Bond modulus of elasticity (N/m^2) 

% Calculate properties of beam cross-section 

area = w.*t; % Beam cross-sectional area (m^2) 

areab — wb*tb; % Bond cross-sectional area (m^2) 

areap = wp*tp; 7o Piezoceramic actuator cross-sectional area (m^2) 


moi = w.*t.^3/12; X *?* Beam cross-section moment of inertia (m^4) 
% Bond cross-section moment of inertia (m4) 
то — wb*tb*(t^2/2 * t*tb -* 2*tb^2/3); 
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% Piezoceramic actuator cross-section moment of inertia (m^4) 

moip — wp*tp*((t*2*tb)^2/2 t (t*2*tb)*tp * 2*tp^2/3); 

EI = E*moitEb*moib+Ep*moip; % Beam stiffness 

K = Ep*d31*wp*(t+2*tb+tp)/EI; % (Moment per unit voltage)/(Beam stiffness) 
K = K*ones(n,1); 


% Set run number, minimum and maximum locations for actuators 
% and spacing between actuator locations in range of positions 
% evaluated between minimum and maximum locations for each actuator 


Run = 6; % Вип number 

ЗауеШе = 'NTII3 6'; 

xmin — [0.055 0.355 0.655]; % Кип 6 (0.005 spacing) 
xmax = [0.095 0.395 0.695]; % Run 6 (0.005 spacing) 


Spacing = [0.005 0.005 0.005]; % Кип 6 


Yo Set initial voltages for unconstrained voltage optimization algorithm 
VO = zeros(n, 1); 


% Set options for unconstrained voltage optimization algorithm 
Options = foptions’; 
Options(14) = 2000; % Maximum number of iterations (0 -> 100*n) 


% Divide each actuator location range into up to 9 equally spaced 
% possible locations separated by the interval set in 'Spacing’ 
Шу == 1 % 1 Actuator 
xrange — [xmin:Spacing:xmax |; 
elseif n — 2 % 2 Actuators 
xrange = [xmin(1):Spacing(1):xmax(1);xmin(2):Spacing(2):xmax(2)]; 
else % 3 Actuators 
xrange = [xmin( 1): Spacing(1):xmax(1);xmin(2):Spacing(2):xmax(2), ... 
xmin(3):Spacing(3):xmax(3)]; 
end 


% Evaluate optimum actuator voltages and beam shape error 
% for each combination of actuator locations 
x = zeros(n,1); 
Maxindex = length(xrange); 
ifn —]1 % 1 Actuator 
Ermatrix = zeros(Maxindex, 1); 
Voltagel = zeros(Maxindex, 1); 
for k = 1:Maxindex 
fprintf([^nk = ' int2str(k) ^n']) 
х = xrange(k); 
V = fmins(‘beamerVu', VO,Options,| ],x,la,K,C,L,n); 
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Voltagel(k) = V; 
Ermatnx(k) = beamerVu(V,x,la,K,C,L,n); 
ifk--1 
if Ermatrix(k) « Minerror 
Optloc = k; 
OptV = V; 
Minerror = Ermatrix(k) 
end 
else 
Optloc = 1; 
OptV = V; 
Minerror = Ermatrix(k) 
end 
end 
clear k 
Optx = xrange(Optloc) 
OptV 
Minerror 


elseif n == 2 Vo 2 Actuators 
Ermatrix — zeros(10); 
Voltagel — zeros(10); 
Voltage2 — zeros(10); 
Voltsum = zeros(10); 
for J = 1:Maxindex 
x(1) = xrange(1,J); 
for k = 1:Maxindex 
fprintf([^njk = ' int2str(]) ' ' int2str(k) ^n']) 
x(2) = xrange(2,k); 
if (x(1)+la(1)<=x(2)) 
V = fmins(‘beamerVu', VO,Options,[],x,la,K,C,L,n); 
Voltagel(J,k) ^ V(1); 
Voltage2(J,k) ^ V(2); 
Ermatrix(j,k) = beamerVu(V,x,la,K,C,L,n); 
if j+k ~= 2 
if Ermatrix(j,k) < Minerror 
Optloc = [j k]; 
OptV = V; 
Minerror = Ermatrix(j,k) 
end 
else 
Optloc - [1 1]; 
OptV = V; 
Minerror = Ermatrix(j,k) 
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end 
else 
Voltagel(j,k) — 0; 
Voltage2(j,k) = 0; 
Ermatrix(j,k) = 0; 
end 
end 
end 
clear j k 
Optx = [xrange(1,Optloc(1));xrange(2, Optloc(2))] 
OptV 
Minerror 
Voltsum = Voltagel+Voltage2; 


else % 3 Actuators 
Ermatrix — zeros(100,10); 
Voltagel — zeros(100,10); 
Voltage2 — zeros(100,10); 
Voltage3 — zeros(100,10); 
Voltsum = zeros(100,10); 
for 1= 1:Maxindex 
x(1) = xrange(1,1); 
for j = 1:Махтдех 
x(2) = xrange(2,j); 
for k = 1:Maxindex 
fprintf([^\nijk = ' int2str(1) ' ' int2str(j) ' ' int2str(k) ^n']) 
x(3) = xrange(3,k); 
if (x(1)+la(1)<=x(2)) & (x(2)+1a(2)<=x(3)) 
V — fmins(beamerVu', VO, Options,[],x,la, K, C,L,n); 
Voltagel(10*(i-1)33,k) ^ V(1); 
Voltage2(10*(i-1)+4,k) = V(2); 
Voltage3(10*(1-1)4+),k) = V(3); 
Ermatrix(10*(i-1)+),k) = beamerVu(V,x,la,K,C,L,n); 
if i+j+k ~= 3 
if Ermatrix(10*(i-1)+),k) < Minerror 
Орџос = [1] k]; 
OptV = У; 
Minerror = Ermatrix(10*(1-1)+),k) 
end 
else 
Optloc = [1 1 1]; 
OptV = V; 
Minerror = Ermatrix(10*(1i-1)+),k) 
end 
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else 
Voltagel(10*(1-1)*j,k) = 0; 
Voltage2(10*(i-1)+j,k) = 0; 
Voltage3(10*(1-1)+),k) = 0; 
Ermatrix(10*(i-1)+j,k) = 0; 
end 
end 
end 
end 
clear 1) k 
Optx = [xrange(1,Optloc(1));xrange(2,Optloc(2));xrange(3,Optloc(3))] 
Opt V 
Minerror 
Voltsum — Voltagel-* Voltage2-* Voltage3; 
end 


% Plot results 
x = linspace(0,L,100); 
y = beamdisp(x,Optx, Opt V,la,n,K); % Fixed length actuators 
yact = beamdisp([Optx;Optx+la],Optx,OptV,la,n,K); 
yf = С*х.^2, 
figure 
plot(x,y,'r-x,yf,'g.' [Optx;iOptx-Ha],yact,'r*^) 
legend('Actual','Desired',- 1) 
grid 
title'Actual vs. Desired Beam Shapes’) 
xlabel(X (m)) 
ylabel('Beam Displacement Y (m)') 
ifn = 1 
text(0.1,0.9,['Run ' int2str(Run) ' Optimum Values for 1 Actuator:'], ... 
"Units ,'Normalized") 
text(0.1,0.8,['x1 =' num2str(Optx(1))'m / V1 =' num2str(OptV())) ... 
' V'|,'Units';,"'Normalized") 
text(0.1,0.7,['Error — ' num2str(Minerror,4) |,'Units','Normalized") 
elseif n — 2 
text(0.1,0.9,["Run ' int2str(Run) ' Optimum Values for 2 Actuators:'], ... 








"Units',"Normalized") 

text(0.1,0.8,['x1 — ' num2str(Optx(1)) ' m / V1 7' num2str(Opt V(1)) ... 
' V'|,'Units','Normalized") 

text(0.1,0.7,[ x2 ^ ' num2str(Optx(2)) ' m / V1 =' num2str(OptV(2)) .. 
' V''Units','Normalized") 


text(0.1,0.6,['Error — ' num2str(Minerror,4) |,'Units',"'Normalized") 
elseif n — 3 
text(0.1,0.9,["Run ' int2str(Run) ' Optimum Values for 3 Actuators:'], ... 
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'Units', Normalized’) 
text(0.1,0.8,[x1 — ' num2str(Optx(1)) ' m / V1 —' num2str(OptV(1)) ... 


' V'|,'Units','Normalized") 

text(0.1,0.7,['x2 = ' num2str(Optx(2)) ' m / V1 =' num2str(OptV(2)) ... 
' V'|,'Units',"Normalized") 

text(0.1,0.6,[ x3 — ' num2str(Optx(3)) ' m / V1 =' num2str(OptV(3)) ... 
' V'|,'Units',"Normalized") 


text(0. 1,0. 5,['Error = ' num2str(Minerror,4)],'Units','Normalized') 
end 


9o Save data to . MAT file 
eval(['save ' Savefile]) 


% PROGRAM TO DETERMINE OPTIMUM ACTUATOR INPUT VOLTAGES 
% FROM A GIVEN DISCRETE SET OF POSSIBLE VOLTAGES 

% 

% This program determines the optimum actuator input voltages 

% from a specified discrete range of possible actuator 

% voltages to achieve the minimum error between the actuator deformed 

% beam shape and a desired parabolic shape function for a cantilever 

% beam with piezoceramic actuators bonded to both sides of the beam. 

% 1, 2 or 3 actuators can be specified 


% Coefficient of desired parabolic shape function 


С = 0.001; 

% Ещег beam properties 

Е = 1; Yo Beam length (m) 

w = 0.0508; % Beam width (m) 

t = 0.001575; % Beam thickness (m) 

Е = 7.2е10, % Beam modulus of elasticity (N/m^2) 
% Enter piezoceramic actuator properties 

% Navy Type II PZT 

n=3; % Number of actuators 

la=0.1*ones(n,1); % Actuator lengths (m) 

wp = 0.0381; % Piezoceramic actuator width (m) 

tp = 0.00027; y/o Piezoceramic actuator thickness (m) 
Ep = 6.1e10; % Piezoceramic actuator modulus of elasticity (N/m^2) 
d31 = 1.71e-10; % Piezoelectric constant (m/V) 

Vmax = 150; % Maximum actuator input voltage (V) 
Vmin = -150; % Minimum actuator input voltage (V) 
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% Enter properties of bond between piezoceramic actuator and beam 


wb = wp; % Bond width (m) (equal to actuator width) 

tb = 0; % Bond thickness (m) 

Eb = 1.78e9; % Bond modulus of elasticity (N/m^2) 

% Calculate properties of beam cross-section 

area = w.*t; 7/o Beam cross-sectional area (m^2) 

areab — wb*tb; ?/o Bond cross-sectional area (m^2) 

areap — wp*tp; Yo Piezoceramic actuator cross-sectional area (m^2) 


moi — w.*t.^3/12; ?/o Beam cross-section moment of inertia (m^4) 

% Bond cross-section moment of inertia (m^4) 

moib — wb*tb*(t^2/2 + t*tb + 2*tb%2/3); 

% Piezoceramic actuator cross-section moment of inertia (m4) 

moip — wp*tp*((t*2*tb)^2/2 + (t+2*tb)*tp + 2*tp%2/3); 

EI = E*mo1t+Eb*moib+Ep*moip; % Beam stiffness 

К = Ep*d31*wp*(t+2*tb+tp)/EI, % (Moment per unit voltage)/(Beam stiffness) 
К = K*ones(n, 1); 


% Set run number, minimum and maximum voltages for actuators 
% and increment between actuator voltages in range of positions 
% evaluated between minimum and maximum locations for each actuator 


Run= 1; Yo Run number 

Savefile ^ 'NTII1 1У"; 

Vmin = [10]; 9o Run 1 (5 V increment) 
Vmax = [50]; % Run 1 (5 V increment) 
Incremnt = [5]; % Кип 1 


% Set actuator locations 


ifn == 1 
x = 0.2420; 
elseif n = 2 
x = [0.1220;0.5394], 
else 
x = [0.0728;0.3678;0.6641]; 
end 


% Set initial voltages for unconstrained voltage optimization algorithm 
VO = zeros(n, 1); 


% Set options for unconstrained voltage optimization algorithm 
Options = foptions’; 


Options(14) = 2000; % Maximum number of iterations (0 -^ 100*n) 


% Divide each actuator voltage range into up to 9 equally spaced 
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% possible locations separated by the interval set in 'Spacing' 
ifn = 1 % 1 Actuator 
Vrange = [ Vmin:Incremnt: Vmax]; 
elseif n — 2 % 2 Actuators 
Vrange — [Vmin(1):Incremnt(1): Vmax(1); Vmin(2):Incremnt(2): Vmax(2)]; 
else % 3 Actuators 
Vrange = 
[ Vmin(1):Incremnt(1): Vmax(1);Vmin(2):Incremnt(2): Vmax(2); Vmin(3):Incremnt(3): Vma 
x(3)]; 


end 


% Evaluate optimum actuator voltages and beam shape error 
Yo for each combination of actuator locations 
Maxindex = length( Vrange); 
if n == 1 % 1 Actuator 
Ermatrix = zeros(Maxindex, 1); 
Voltage] = zeros(Maxindex, 1); 
for k = 1:Maxindex 
fprintf([^nk — ' int2str(k) ^n']) 
V — Vrange(k); 
Ermatrix(k) = beamerVu(V,x,la,K,C,L,n); 
ifk ~= 1 
if Ermatrix(k) < Minerror 
OptV = V; 
Minerror = Ermatrix(k) 
end 
else 
OptV = V; 
Minerror = Ermatrix(k) 
end 
end 
clear k 
OptV 
Minerror 


elseif n = 2 % 2 Actuators 
Ermatrix — zeros(10); 
for J » 1:Maxindex 
for k = 1:Maxindex 
fprintf([Anjk — ' int2str(j) '' int2str(k) ^n']) 
V — [Vrange(1,j); Vrange(2,k)]; 
Ermatrix(j,k) = beamerVu(V,x,la,K,C,L,n); 
if jik ~= 2 
if Ermatrix(j,k) « Minerror 
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OptV = У; 
Minerror = Ermatrix(j,k) 
end 
else 
OptV = V; 
Minerror — Ermatrix(j,k) 
end 
end 
end 
clear J k 
OptV 
Minerror 


else % 3 Actuators 
Ermatrix — zeros(100,10), 
for 1— 1: Maxindex 
Гог ] = 1:Махшдех 
for k = 1:Maxindex 

fprintf([^nijk — ' int2str(1) '' int2str(j) '' int2str(k) ^n']) 
V — [Vrange(1,1); Vrange(2,]); Vrange(3,k)]; 
Ermatrix(10*(i-1),k) = beamerVu(V,x,la,K,C,L,n); 


if Hk — 3 
if Ermatrix(10*(1-1)43,k) < Мшеггог 
OptV = V; 


Optindex = [i J К]; 
Minerror — Ermatrix(10*(1-1),k) 
end 
else 
OptV = У; 
Optindex - [1j k]; 
Minerror = Ermatrix(] 0*(1-1)+),k) 
end 
end 
end 
end 
clear ijk 
OptV 
Minerror 
end 


Yo Plot results 

ifn=1 %1 Actuator 
% Plot variation in shape error vs. actuator voltage 
figure 
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plot(Vrange,Ermatrix) 

title(Error vs. Actuator Input Voltage for 1 Actuator’) 

xlabel(‘Actuator Voltage (V)') 

ylabel(‘Error') 

text(0.15,0.95,['Run ' int2str(Run)],'Units',"'Normalized") 
text(0.15,0.85,|'Actuator Location: ' num2str(x) ' m'],'Units',"'Normalized") 
grid 


elseifn — 2 %2 Actuators 
Yo Plot variation in shape error vs. actuator positions 
figure 
meshz(Vrange(2,:), Vrange( 1 ,:),Ermatrix( 1: Maxindex, ]: Maxindex)) 
title(‘Error vs. Actuator Input Voltages for 2 Actuators’) 
xlabel( Voltage 2 (V)') 
ylabel('Voltage 1 (V)') 
text(0.15,0.95,['Run ' int2str(Run)],"Units',"Normalized’) 
text(0.15,0.85,['Actuator Locations: 'num2str(x(1))'m '... 
num2str(x(2)) ' m'],"Units','Normalized") 


else % 3 Actuators 
% Plot variation in shape error for Actuator 1 fixed 
% at optimum location 
figure 
meshz( Vrange(3,:), Vrange(2,:),Ermatrix(10*(Optindex(1)-1)+1: ... 
10*(Optindex(] )-1)+Maxindex, | :Maxindex)) 
title(Error vs. Actuator 2 & 3 Voltage at Optimum Actuator 1 Voltage’) 
xlabel('Voltage 3 (V)') 
ylabel('Voltage 2 (V)') 
text(0.15,0.95,['Run ' int2str(Run)],"Units',"Normalized’) 
text(0.15,0.85,['Actuator Locations: ' num2str(x(1)) ' m ' ... 
num2str(x(2)) ' m ' num2str(x(3)) ' m'],"Units',"Normalized") 


% Plot variation in shape error for Actuator 2 fixed 

% at optimum location 

figure 

meshz(Vrange(3,:), Vrange(1,:),Ermatrix(] 0*[0:Maxindex-1]+ ... 
Optindex(2)*ones(1,Maxindex), 1: Maxindex)) 

title(‘Error vs. Actuator 1 & 3 Voltage at Optimum Actuator 2 Voltage’) 

xlabel("Voltage 3 (V)') 

ylabel('Voltage 1 (V)') 

text(0.15,0.95,['Run ' int2str(Run)],'Units','Normalized') 

text(0.15,0.85,|'Actuator Locations: ' num2str(x(1)) ' m ' ... 
num2str(x(2)) ' m ' num2str(x(3)) ' m'],'Units',"'Normalized') 
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y/o Plot variation in shape error for Actuator 3 fixed 

% at optimum location 

figure 

Ermesh = reshape(Ermatnx(:,Optindex(3)), 10,10)’; 

meshz(Vrange(2,:), Vrange(1,:),Ermesh(] :Maxindex, 1: Maxindex)) 

title(Error vs. Actuator 1 & 2 Voltage at Optimum Actuator 3 Voltage’) 

xlabel('Voltage 2 (V)') 

ylabel('Voltage 1 (V)? 

text(0.15,0.95,['Run ' int2str(Run)],'Units','Normalized") 

text(0.15,0.85,['Actuator Locations: ' num2str(x(1)) ' m ' ... 

num2str(x(2)) ' m ' num2str(x(3)) ' m'],'Units;'Normalized") 

end 


% Save data to .MAT file 
eval(['save ' Savefile]) 
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